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Abstract 


The  mathematical  modeling  of  electrical  and  thermal 
phenomena  that  occur  during  in  situ,  low  frequency, 
electrical  heating  of  oil  sand,  or  other  materials  with 
similar  electrical  properties,  is  studied.  Particular 
attention  was  directed  toward  the  modeling  of  the  electrical 
field  at  an  interface  between  two  media  of  different 
electrical  conductivities.  The  finite  difference  method  is 
used  to  develop  a  computer  program  to  simulate,  in  two 
dimensions,  the  mathematical  model  that  was  developed. 

This  computer  program  is  used  to  study  the  role  of 
various  parameters  in  determining  the  final  temperature 
profile  after  a  fixed  time  period  and  fixed  rate  of  heating. 
The  parameters  studied  include  the  position  of  the 
electrodes  relative  to  the  overlying  and  underlying 
formations,  the  spacing  between  the  electrodes,  and  the 
ratio  of  the  conductivity  of  the  oil  sand  formation  to  the 
conductivity  of  the  surrounding  formations. 
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1.  Introduction 


1.1  Electrical  Preheat  In  Situ  Method 

In  northern  Alberta  there  exist  large  deposits  of 
petroleum  in  the  form  of  oil  sand.1  The  oil  in  these 
deposits  is  too  viscous  to  recover  by  conventional  methods. 
Commercial  developments  to  date  have  mined  the  oil  sand  in 
open  pit  mines,  then  separated  out  the  petroleum  (in  the 
form  of  bitumen)  in  large  extraction  plants.  The  bitumen  is 
upgraded  to  "synthetic  crude",  then  piped  to  a  refinery. 

As  the  thickness  of  the  earth,  or  overburden,  above  the 
oil  sand  formation  increases,  open  pit  mining  becomes 
uneconomical.  Only  about  ten  percent  of  the  oil  sand 
reserves  are  recoverable  by  surface  mining.2  Thus,  there  is 
considerable  interest  in  developing  in  situ  (in  place) 
methods  of  separating  out  the  bitumen.  A  wide  range  of 
possible  in  situ  techniques  have  been  suggested.3 

One  of  the  proposed  in  situ  techniques  is  steam 
flooding.  In  the  steam  flood  method  an  array  of  injection 
and  production  wells  are  drilled,  then  steam  is  pumped  down 
the  injection  wells  and  oil  and  water  are  pumped  up  the 
production  wells.  A  major  problem  with  this  method  is  in 
developing  "communication"  between  the  injection  and  the 


production 

wells. 

The 

unheated 

bi tumen 

is  so  vi scous  i t 

blocks  the 

steam 

and 

heated , 

mob i 1 i zed 

bitumen  from 

traveling  from  one  well  to  the  other.  Various  methods  have 
been  proposed  to  develop  communication.  These  include 


1 
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injection  of  solvents,  formation  of  emulsions  and  formation 
f ractur i ng . 1 7  The  first  two  methods  are  designed  to  reduce 
the  viscosity  of  the  bitumen,  and  the  third  method  produces 
a  direct  channel  between  wells. 

The  electrical  preheat  in  situ  method  proposes  that  the 
oil  sand  formation  be  heated  by  passing  an  electrical 
current  through  it,  and,  once  the  vicosity  has  been  lowered 
sufficently,  steam  or  some  other  hot  fluid  be  used  to  drive 
out  the  bitumen.  The  viscosity  of  bitumen  decreases  rapidly 
with  increases  in  temperature.4  At  natural  formation 
temperatures  Athabasca  bitumen  has  a  viscosity  in  the 
millions  of  centipoise.  If  this  bitumen  is  heated  to  80*C 
the  viscosity  drops  to  1000  cp,  and  if  heated  to  200* C  the 
viscosity  will  be  about  10  cp.  The  heat  used  to  raise  the 
temperature  of  the  bitumen  is  generated  right  in  the  oil 
sand  by  the  electrical  current  flowing  in  the  connate  water. 
This  heating  in  place  overcomes  some  of  the  problems  usually 
associated  with  heating  heavy  oil  formations  with  steam, 
such  as  low  thermal  conductivity  and  difficulty  directing 
the  applied  heat.  Electrical  heating  has  seen  limited  use  in 
secondary  recovery  of  conventional  oil  fields.5  Several 
patents  have  been  issued  on  the  application  of  electrical 
heating  to  oil  sand.  Petro-Canada  Explorations  Ltd.  is 
running  a  pilot  plant  investigating  the  electrical  preheat 
method . 4 


% 


3 


1.2  Means  of  Investigation 

Before  pilot  plant  tests  or  commercial  developments  of 
the  electrical  preheat  method  are  embarked  upon  it  is 
necessary  to  develop  an  idea  of  the  optimal  electrode 
configuration,  heating  rates,  and  heating  times.  Other  than 
pilot  plant  or  field  tests,  the  study  of  the  electric 
preheat  method  may  be  carried  out  by  three  means:  analytic 
calculations,  laboratory  scale  models,  and  numerical 
simulation.  In  all  three  of  these  means  it  is  necessary  that 
a  mathematical  model  of  the  important  processes  be 
deve 1 oped . 

The  electrical  conductivity  of  oil  sand  is  temperature 
dependent.  This  couples  the  governing  thermal  and  electrical 
equations.  The  resulting  set  of  partial  differential 
equations  is  difficult  or  impossible  to  solve  analytically 
for  all  but  the  simplest  one  dimensional  cases. 

Physical  scale  modeling  of  in  situ  recovery  methods  is 
discussed  by  Farouq  Ali  and  Redford,6  and  the  physical 
modeling  of  electromagnetic  heating  of  oil  sand  is  discussed 
by  Vermeulen,  Chute,  and  Cervenan.7 

It  was  decided  to  use  numerical  simulation  to  study 
some  aspects  of  the  electrical  preheat  method.  Two  other 
studies  reporting  numerical  simulation  of  electrical  heating 
of  petroleum  reservoirs  are  by  Todd  and  Howell8  and 
El-Feky.9  Todd  and  Howell  simulate  electrically  heating  oil 
sand  to  several  hundred  degrees  Celsius,  while  cooling  the 
well  bore.  El-Feky  studied  electrically  heating  conventional 


4 


oil  reservoirs  for  secondary  recovery.  Both  of  these  studies 
and  simulators  differ  in  detail  and  application  from  this 
study. 

Numerical  simulation  may  be  used  as  a  compliment  or 

alternative  to  physical  scale  modeling.  Numerical  simulation 

allows  easy  adjustment  of  parameters  and  electrode 

configurations,  but  is  only  as  good  as  the  mathematical 

model  on  which  the  simulation  is  based.  Comparison  of 

numerical  simulation  results  and  physical  scale  model 

results  can  provide  insight  into  the  accuracy  of  a 

mathematical  model.  Physical  scale  modeling  of  the 

electrical  heating  of  oil  sand  is  being  done  by  Vermeulen, 

Chute  and  Cervenan7  at  the  University  of  Alberta,  and  some 

of  the  results  of  their  group  are  compared  to  the  results 

obtained  by  the  numerical  simulator  which  is  discussed  in 

# 

this  thesis. 


1.3  Electrical  and  Thermal  Properties  of  Oil  Sand 

The  electrical  conductivity  and  relative  dielectric 
constant  of  reconstituted  oil  sand  were  measured  by  Chute, 
Vermeulen,  Cervenan  and  McVea10  and  correlated  to  frequency, 
moisture  content,  temperature,  and  density.  Similarly,  the 
heat  capacity  and  thermal  conductivity  of  various  grades  of 
oil  sand  were  reported  by  Cervenan,  Vermeulen,  and  Chute11. 
The  applicable  results  of  these  two  papers  are  summarized  in 
Table  1.1  and  Table  1.2. 


. 
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Table  1.1  Electrical  Properties  of  Reconstituted  Oil  Sand  at 

60  Hz  and  24*C. 10 


Material  Water 

(%  wt. ) 

Oi  1  Sand  A  5.8 

Oi 1  Sand  B  3.3 

Oi 1  Sand  C  1.3 


Conductivity  Dielectric 

(S/m)  Constant 


1.5X10-2  4.X  104 
6.  X  10-3  1 .  X  104 
1.2  X  IQ'3  1 .  X  103 


Table  1.2  Thermal  Properties  of  Reconstituted  Oil  Sand  and 

Other  Materials.11  18 


Material 

Water 

Conduct ivi ty 

Heat  Capacity 

(%  wt. ) 

(W/m-K) 

( d/m^-K ) 

Oi 1  Sand 

4.3 

1.80 

Oi 1  Sand 

2.1 

1.42 

- 

Oi 1  Sand 

1.3 

1.35 

- 

Oi 1  Sand 

1.5 

- 

1.71  X  106 

Oi 1  Sand 

11.1 

- 

2.23  X  106 

Oil  Sand 

1.4 

— 

1 .83  X  106 

Sandstone 

- 

.877 

1.59  X  106 

Limestone 

- 

1.70 

1.86  X  106 

Shale 

- 

1 .04 

1.87  X  106 
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The  electrical  conductivity  was  found  to  increase 
approximately  linearly  with  temperature  dependence  up  to 
90  C.  The  computer  program  written  in  this  study 
incorporates  this  temperature  dependence.  The  electrical 
conductivity  is  expressed  by 

O'(T)=0J,  [l+oc(T-24)]  (1.1) 

where  T  is  in  degrees  Celsius  and  C/  is  in  Siemens/meter . 
The  constant  oc  has  units  *C_1  and  is  the  temperature 
dependence  of  electrical  conductivity.  A  typical  measured 
value  of  oC  for  oil  sand  is  2.3  X  10"2  ’C"1. 

The  author  knows  of  no  published  data  on  the 
temperature  dependence  of  conductivity  of  oil  sand  at 
temperatures  above  90* C.  Chute  et.al.10  report  that  above 
90 ‘C  the  conductivity  tends  to  level  off,  and  hence  the 
value  of  oc  decreases.  In  an  electrical  preheat  of  an  oil 
sand  formation  the  temperature  will  remain  below  90*C  in  the 
bulk  of  the  formation.  Near  the  electrodes  the  temperature 
may  reach  250* C  or  higher,  but  the  conductivity  of  the  oil 
sand  will  probably  be  affected  by  the  injection  of  brine 
into  the  formation.  Due  to  the  lack  of  published  data  the 
author  has  assumed  that  the  conductivity  increases  linearly 
with  temperature,  even  well  above  90* C. 

While  the  thermal  conductivity  and  heat  capacity  of  oil 
sand  may  have  some  temperature  dependence,  there  is  no 
detailed  study  on  such  temperature  dependence  reported  in 
the  literature.  A  few  measurements  by  Cervenan  et.  al . 1 1 


- 
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show  no  marked  temperature  dependence  for  the  heat  capacity 
of  oil  sand  samples  between  20°C  and  70#C.  The  computer 
program  was  written  so  that  the  temperature  dependences  of 
thermal  conductivity  and  heat  capacity  could  easily  be 
included  once  these  dependences  are  known. 

The  effects  on  the  electrical  preheat  process  of 
changing  the  ratio  of  the  electrical  conductivities  of  the 
oil  sand  and  surrounding  formations  is  one  of  the  items 
studied  in  this  thesis.  For  the  study  of  the  effects  of 
changing  this  ratio  the  thermal  properties  of  the 
surrounding  formations  are  assumed  to  be  the  same  as  those 
of  the  oil  sand.  It  is  also  possible  to  give  the  surrounding 
formations  different  thermal  properties  than  the  oil  sand. 


1.4  The  Computer  Program  "MEGAERA" 

A  computer  program  entitled  MEGAERA16  was  written  to 
numerically  simulate  the  electrical  conduction  heating  of 
.oil  sands,  or  other  solid  materials.  The  mathematical  model 
and  finite-  difference  scheme  used  in  MEGAERA  are  discussed 
in  detail  in  the  second  and  third  chapters  of  this  thesis. 
MEGAERA  was  used  to  study  a  variety  of  electrical  heating 
configurations  with  different  conductivity  ratios  between 

the  oil  sand  and  the  surrounding  formations.  The  results  of 

« 

this  study  are  given  in  chapter  four  of  this  thesis. 

Any  study,  analytic,  experimental  or  numerical,  must  be 
limited  to  some  extent.  The  specification  for  MEGAERA  were 
based  on  the  type  of  physical  scale  modeling  being  done  by 
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Vermeu len ,  Chute  and  Cervenan  at  the  University  of  Alberta 
and  the  experience  of  the  author  in  writing  and  using  an 
earlier  two  dimensional  electrical  heating  simulator. 

There  were  two  primary  considerations  in  developing 
MEGAERA.  First,  the  program  had  to  be  able  to  simulate 
overburden  and  underburden  of  different  electrical 
conductivities  than  the  oil  sand.  The  optimal  electrode 
configuration  for  heating  an  oil  sand  formation  is  highly 
dependent  on  whether  the  surrounding  formations  are  more  or 
less  conductive  than  the  oil  sand.  Second,  the  program  must 
be  relativly  inexpensive  to  run.  For  this  reason  MEGAERA  was 
limited  to  two  dimensions.  A  typical  MEGAERA  run  takes  two 
to  three  minutes  of  C.P.U.  time  on  an  Amdahl  470V/7.  The 
time  required  by  a  three  dimensional  program  can  be 
conservat i vly  estimated  by  multiplying  the  time  required  by 
the  two  dimensional  program  by  the  number  of  grid  lines  in 
the  third  dimension.  For  a  reasonable  amount  of  detail  in 
the  third  dimension  the  program  would  be  too  expensive  to 
run  more  than  a  few  times. 

While  the  universe  is  generally  perceived  to  have  three 
spacial  dimensions,  it  is  often  convenient  and  reasonably 
accurate  to  model  a  physical  system  as  two  dimensional .  In 
studying  electrical  heating  of  oil  sand  there  are  two 
different  two  dimensional  slices  of  interest.  A  horizontal 
(areal)  study  of  the  heating  pattern  in  a  five  spot  pattern 
assumes  that  the  oil  sand  formation  is  infinitly  thick,  but 
may  yield  results  which  are  correct  for  the  middle  of  a 
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the  simple  cylindrical  variety, 
would  be  used . 


Figure  1 . 1 
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finite  formation.  A  vertical  section  studying  the  heating 
between  two  parallel  plates  allows  one  to  see  the  effects  of 
changing  the  conductivity  of  the  surrounding  formations. 
Parallel  plates  may  be  approximated  in  a  field  test  by 
systems  of  closely  spaced  electrodes  or  by  using  a  more 
widely  spaced  system  of  electrodes  in  which  each  electrode 
has  been  effectively  enlarged  by  jet  cutting  the  formation 
with  a  brine  solution. 

In  modeling  configurations  with  overburden  and 
underburden  the  energy  losses  due  to  electrical  current 
leakage  and  heat  conduction  from  the  oil  bearing  formation 
are  calculated  by  MEGAERA.  Voltage,  current,  power  and 
resistance  calculation  are  done  at  each  time  step.  MEGAERA 
will  automatically  scale  the  applied  voltages  to  produce 
constant  power  and  constant  current  heating,  if  this  is 
desired.  A  post  processor,  called  MEGAERA . PLOTTER ,  was 
written  to  produce  contour  and  profile  plots  of  the 
electrical  potential,  heating  rate  and  temperature  at 
specified  times. 


2.  The  Mathematical  Model 


2.1  Electrical  Fields 

The  basic  set  of  equations  of  classical  electricity  and 
magnetism,  first  written  by  Maxwell  in  1863,  are:12 


^•D7=  (2.D 

7-6=0  (2.2) 

(2-3) 

(2.4) 


These  equations,  together  with  the  appropriate  boundary 
conditions,  may  be  used  to  find  the  electric  field  of  any 
electrical  heating  configuration  at  any  frequency.  However, 
an  analytic  or  numerical  solution  to  Maxwell's  equations  is 
difficult,  if  not  impossible,  in  all  but  the  simplest  cases. 
It  is  therefore  desirable  to  simplify  these  equations, 
dropping  terms  which  are  negligible  for  the  frequencies  and 
electrode  configurations  of  interest. 

2.1.1  The  Quasi -static  Approximat ion 

In  their  paper  on  physical  scale  modeling,  Vermeulen 
et.al.7  discuss  the  simplification  of  Maxwell's  equations 
for  various  types  of  electrical  heating.  The  discussion  in 
this  section  (2.1.1)  is  drawn  from  the  work  of  their  paper 


12 


and  from  Magid . 1 3 

Consider  the  displacement  current  term  in  equation 
(2.4).  The  ratio  of  conduction  current  to  displacement 
current  is  given  by  ^AjjS  for  time  harmonic  solutions.  For 
oil  sand  this  ratio  is  typically  greater  then  one  hundred  at 
power  frequencies  (60  Hz),  and  does  not  fall  to  ten  until 
frequencies  of  about  one  megahertz  are  reached.  Thus,  for 
power  frequencies  the  effects  of  displacement  current  may  be 
neglected. 

The  study  of  the  fields  of  two  dimensional,  parallel 
electrode,  electrical  heating  configurations  is  similar  to 
the  calculation  of  electric  fields  in  lossy  parallel  plate 
capacitors.  In  these  types  of  problems  the  effects  of  the 
time  varying  magnetic  field  on  the  time  varying  electric 
field  may  be  neglected  if  the  largest  dimension  of  the 
problem  geometry  is  much  less  than  a  wavelength  in  the 
medium.  This  reduction  of  Maxwell's  equations  to  the  static 
field  equations  by  neglecting  these  magnetic  and 
displacement  current  effects  is  known  as  the  quasi -static 
approximation. 

The  wavelength  in  oil  sand  at  60  Hz  is  between  2000  m 
and  9000  m,7  depending  on  the  water  content  of  the  oil  sand. 
As  the  maximum  dimension  of  a  typical  electrical  heating 
configuration  will  probably  be  less  then  200  m,  the 
quasi -static  approximation  should  be  valid.  If  the 
surrounding  formations  (overburden  and  underburden)  are  also 
to  be  modeled,  the  wavelength  in  these  materials  must  also 
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be  much  greater  then  the  maximum  problem  dimension  and  the 
loss  tangent  must  be  much  greater  than  one. 

When  using  the  quasi -static  approximation  the  geometry 
of  the  electric  field  is  calculated  from  the  static  electric 
field  equations.  While  calculating  the  fields  the  potential 
difference  between  electrodes  is  fixed  at  the  peak  A.C. 
value.  The  time  varying  electric  field  is  then  expressed  as: 

t  -  E(x,y,z)  5 in  uut  (2.5) 

where  El(xt/,z)  describes  the  static  field,  u/-2irf  where  f 
is  the  frequency,  and  t  is  the  time  in  seconds.  The  r.m.s. 
value  of  the  potential  difference  is  used  in  MEGAERA  rather 
then  the  peak  value.  This  simplifies  the  calculation  of 
heating  rates  and  conforms  to  the  practise  of  stating  power 
line  voltages  as  r.m.s.  rather  then  peak  values. 

The  discussion  in  the  previous  paragraph  applies  to 
single  phase  heating.  When  three  phase  heating  is  used  the 
geometry  of  the  electric  field  of  each  phase  must  be 
calculated  separately,  using  the  static  field  equations  and 
setting  the  electrodes  .of  the  other  phases  to  ground 
potential.  Each  of  the  three  electric  field  geometries  thus 
calculated  is  then  multiplied  by  the  appropriate  time 
variation,  which  must  include  the  phase  differences. 
Finally,  the  three  fields  are  added  vectorial ly  and  the 
square  of  the  magnitude  of  the  resulting  field  is  averaged 
over  a  cycle  to  get  the  heating  rate.  This  procedure  was  not 
incorporated  in  MEGAERA,  so  multiphase  heating  cannot  be 
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simulated. 

2.1.2  Current  Continuity  Equation 

Neglecting  the  effect  of  the  time  varying  magnetic 
field  of  the  electric  field  yields: 

(2.6) 

Thus,  a  scalar  potential  ^(x,y,z)  may  be  defined.14 

E*  -grad  ^  (2.7) 

Droping  the  displacement  current  term  in  equation  (2.4) 
gives 

V  *  J  (2.8) 

Taking  the  divergence  of  both  sides  of  equation  (2.8)  gives 

divJ-0  (2.9) 

which  states  that  the  current  is  continuous.  As  T  -  crE, 
where  o'  is  the  electrical  conductivity,  equation  (2.9) 
becomes 

div(cr  jrod  -0  (2.10) 

As  the  electrical  conductivity  O'  depends  on  temperature, 
it  will,  in  general,  be  a  function  of  of  position,  and 
equation  (2.10)  may  not  be  simplified  by  moving  the 
conductivity  term  outside  the  divergence  operator. 
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2.1.3  Electrical  Boundary  Conditions 

Three  different  types  of  electrical  boundaries  are  to 
be  modeled  in  the  electrical  conduction  heating  simulation. 
They  are:  electrode  -  conducting  material  boundaries, 
insulator  -  conducting  material  boundaries,  and  boundaries 
between  two  conducting  materials. 

Electrodes  are  simply  modeled  as  regions  of  constant 
voltage.  If  constant  power  or  constant  current  heating  is 
requested,  MEGAERA  first  calculates  the  fields  assuming 
constant  voltage,  then  measures  the  current  and  scales  the 
electric  field  to  the  appropriate  case.  Electrode  -  Oil  Sand 
contact  impedances  were  noted  by  Chute  et.  al.10  when 
measuring  the  electrical  properties  of  oil  sand.  A  major 
cause  of  contact  impedence  is  poor  physical  contact  between 
the  oil  sand  and  the  electrode.  In  the  mathematical  model 
used  it  is  assumed  that  contact  impedance  has  been  made 
negligible  by  a  suitable  technique  such  as  the  injection  of 
a  saline  solution  near  the  electrode. 

Insulator  -  conducting  material  boundaries  are  modeled 
in  electrostatics  as  surfaces  at  which  there  is  no  current 
flow  normal  to  the  surface  (  a  homogeneous  Neumann  boundary 
condition).  Having  no  current  normal  to  the  surface  implies 
that  the  normal  component  of  the  electric  field  is  zero  just 
inside  the  conducting  material.  As  well  as  modeling 
insulating  boundaries,  this  type  of  boundary  condition  may 
be  used  to  model  planes  of  symmetry.  Many  of  the  electrode 
configurations  of  interest  in  electrical  conduction  heating 
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of  oil  sand  consist  of  large  arrays  of  electrodes  of 
alternating  polarity.  The  planes  of  symmetry  in  the  electric 
field  may  be  used  to  divide  the  large  array  into  many  small 
identical  cells.  The  simulation  of  the  heating  in  one  of 
these  small  cells  is  all  that  is  necessary  to  find  the 
heating  in  the  entire  array. 

At  a  boundary  between  two  media  of  finite  conductivity 
the  electric  field  boundary  conditions  are:15 


n  *  (E2-  Z,)  =  0 

(2.11) 

n-(o2-0, )  =  P% 

(2. 12) 

where  E,  ,  Ej  ,  and  D,  ,  D*  ,  are  the  fields  in  the  two  media, 
and  where  is  the  surface  charge  at  the  interface  in 
coulombs  per  square  meter  and  n  is  a  unit  vector  locally 
normal  to  the  boundary.  These  are  general  boundary 
conditions  for  Maxwell's  equations  and  are  valid  for  all 
frequencies.  The  boundary  conditions  may  be  expressed  in 
terms  of  the  normal  and  tangential  components  of  the 
electric  field. 

Starting  with  Maxwell's  equations,  take  the  divergence 
of  equation  (2.4)  and  substitiute  in  equation  (2.1)  to  get 
the  current  continuity  equation. 

V-?  +§ £  =0  (2.13) 

Writing  this  equation  in  its  integral  form  and  applying  the 
divergence  theorem  gives 
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Figure  2.1  Planes  of  Symmetry  in  a  Five  Spot  Array 
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J[  ds  =  ~l  It  dv 


(2. 14) 


To  apply  equation  (2.14)  to  the  boundary  between  two 
materials  of  finite  conductivity,  consider  a  small 
cyl indr i cal ly  shaped  volume  enclosing  a  portion  of  the 
boundary.  One  face  of  the  cylinder,  denoted  St  ,  is  in 
material  1  and  is  parallel  to  the  boundary.  The  second  face, 
S2  ,  is  in  material  2  and  is  parallel  to  the  boundary.  The 
side  of  the  cylinder,  S3  ,  cuts  across  the  boundary 
connecting  the  two  faces.  Applying  equation  (2.14)  now  gives 


(2.15) 


Now,  let  the  length  of  the  side  of  the  cylinder  shrink 
toward  zero.  As  both  materials  1  and  2  are  of  finite 
electrical  conductivity,  the  current  density  must  remain 
finite,  and  the  integral  over  ^  goes  to  zero.  For  S,  and  S2 
sufficently  small  (but  an  order  larger  then  S3  ),  equation 

( 2 . 15)  becomes 


(2.16) 


where  ^  S  is  the  area  of  a  face  of  the  cylinder,  is  the 


surface  charge  density  per  unit  area,  and  £n/  and  Eng  are 


the  components  of  the  electric  field  normal  to  the  boundary 
in  materials  1  and  2  respectively.  Dividing  out  the  *S  and 
assuming  a  time  harmonic  field  yields 


or  E„,  -  05  Zni  =  -j  <Y>S 


(2. 17) 


From  equation  (2.12)  and  assuming  that  the  electrical 
permittivity  is  uniform  and  isotropic  in  each  material 


(2.18) 


Substituting  in  equation  (2.17)  gives 


(o',  )  Eni  -  (cl  +jcu£z)  En 2 


(2.19) 


Recall  that  for  the  quasistatic  approximation  to  hold  G'Aui 
must  be  much  greater  then  one.  Thus,  when  using  the 
quasi  static  approximation  the  normal  component  of  current  is 
continuous  across  a  boundary  between  two  materials  of  finite 
electrical  conductivity. 

Treating  the  problem  as  quasistatic  and  assuming  there 
are  no  double  layers  of  charge  at  the  interface,  the 
electric  potential  must  be  continuous  across  the  boundary. 
Thus  the  boundary  conditions  used  are 


Em  “  OJ  Ena 


(2.20) 


(2.21) 


In  this  discussion  of  the  boundary  between  two 
materials  of  finite  conductivity  it  was  assumed  that  the 
conductivity  was  discontinuous  at  the  boundary.  In  an  actual 
oil  sand  formation  the  conductivity  may  change  continuously 
between  the  oil  sand  and  the  overburden  or  underburden.  The 
electric  fields  calculated  assuming  either  a  sharp, 


discontinuous  boundary,  or  assuming  a  continuous  transition 
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over  a  short  distance  will  be  the  same  except  in  the 
transition  region.  As  no  data  is  available  describing  how 
the  conductivity  changes  from  the  oil  sand  to  the 
surrounding  formations,  and  as  the  transition  distance  is 
likely  small  compared  to  the  formation  thickness,  it  is 
assumed  that  a  sharp,  discontinuous  change  in  conductivity 
occurs  at  the  boundary. 


2.2  Heat  Generation  and  Tranfer 

The  instantaneous  heating  rate  due  to  electrical 
conduction  is  given  by 

<?  =  O'  i-E  (2.22) 

If,  when  using  the  quasi -static  approximation,  the  r.m.s. 
values  of  the  applied  voltages  are  used,  then  the  average 
heating  rate  is 

Qav  =  cr|Ermsf  (2.23) 

Of  the  three  major  heat  transfer  mechanisms  (radiation, 
convection,  and  conduction)  only  conduction  is  included  in 
the  mathematical  model.  Radiation  is  negligible  as  the 
radiation  from  a  black  body  at  reservoir  temperatures  would 
only  penetrate  a  very  short  distance  into  the  oil  sand. 
Natural  convection  due  to  temperature  gradiants  is  neglected 
even  though  its  importance  as  a  heat  transfer  mechanism  in 
oil  sand  at  the  temperatures  of  interest  is  unclear.  The 
scale  modeling  factors  developed  by  Vermeulen  et.  al.7 
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included  only  conduction  in  their  derivation.  By  comparing 
the  results  of  numerical  simulations  based  on  this  premise 
with  the  results  of  physical  scale  model  runs  some  grasp  on 
the  accuracy  of  this  mathematical  model  may  be  obtained.  The 
treatment  of  forced  convection  due  to  fluid  injection  or 
pumping  is  beyond  the  scope  of  this  study. 

Thus,  the  equation  used  to  model  the  heat  generation 
and  transfer  is 

/VlJJ  =  V-fa  VTj  +c'lErJ  (2.24) 

where  /W  is  the  volumetric  heat  capacity,  khis  the  thermal 
conductivity,  and  T  is  the  temperature  in  Celsius.  Note  that 
in  the  computer  program  the  volumetric  heat  capacity  and 
thermal  conductivity  are  assumed  constant  within  each 
formation . 

The  outer  boundaries  are  assumed  to  be  insulated, 
no-heat-flow  boundaries.  If  the  boundaries  are  planes  of 
symmetry  or  actual  insulated  boundaries  this  is  a  good 
model.  However,  if  the  problem  to  be  modeled  consists  of  an 
oil  sand  formation  surrounded  above  and  below  by  formations 
which  are  much  thicker  then  the  oil  sand  formation  (i.e. 
semi - i nf i ni te  overburden  and  underburden)  artificial 
boundaries  will  have  to  be  introduced  to  limit  the  extent  of 
the  problem  domain.  If  these  boundaries  are  far  enough  from 
the  electrodes  little  current  or  heat  will  reach  the  area 
near  these  artificial  boundaries. 

Different  materials  in  the  problem  domain  may  have 
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different  thermal  properties.  At  the  boundary  between  two 
different  materials  the  temperature  and  heat  flow  across  the 
boundary  are  assumed  to  be  continuous. 


3.  Numerical  Solution  of  the  Model 
Finite  difference  techniques  are  used  to  approximate  the 
solution  of  partial  differencial  equations  in  many  different 
fields,  including  oil  reservoir  simulation,  laser-plasma 
interaction  studies,  and  water  flow  problems.  Finite 
difference  approximations  to  derivatives  were  used  by  L. 
Euler  (1768)  and  there  are  currently  several  textbooks 
available  which  explain  the  finite  difference  method  in 
detail.19  20  21  The  error,  stability  and  convergence  of  the 
approximation  and  solution  techniques  are  areas  of  concern 
and  study  when  using  this  method  t‘o  solve  partial 
differential  equations.  In  this  chapter  the  particular 
methods  used  in  developing  MEGAERA  are  outlined,  and  an 
explanation  of  how  boundary  conditions,  stability,  and 
convergence  problems  were  delt  with  is  given. 


3.1  Differencing  of  the  Current  Continuity  Equation 

The  electrical  current  continuity  equation  (equation 
2.10)  may  be  rewritten  (for  two  dimensions)  as 


The  finite  difference  approximation  for  the  first  term  of 
equation  (3.1)  is  given  by 


(3.2) 
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See  figure  3.1  for  details  on  the  grid  structuring  and 
1 abe ling. 

This  approximation  (equation  3.2)  is  "conservative"  in 
that  the  current  flowing  out  of  the  left  boundary  of  the 
grid  block  (i,j)  is  equal  to  the  current  flowing  into  the 


right  boundary  of  the  grid  block  (i+1,j).  A  symmetric 

approximation  for  the  second  (y  direction)  term  of  equation 

(3.1)  is  added  to  equation  (3.2)  to  give  an  algebraic 

difference  equation  at  the  center  of  each  grid  block. 

Rearranging  the  terms  gives 

R’j  *  c'j  t-'-j  *  Fu  A h  * 0 

(3.3) 

where 

R/j  +  or/.j  )/l  kj  ( ks  *  kj.J] 

0 

(3.4) 

S,j  1  ( <T;„j  *  o;-'  )/C  h;  ( h;  +  h;.,  J] 

(3.5) 

A"  -  -  —  S • .  —  C'*  ~  P' * 

(3.6) 

C,j-  •(er,-.,ij  +  <X-j  )/i  h,-.,  (h;  *h;., 

(3.7) 

&j)/[  kj.,  (  kj  *  kj„ )] 


(3.8) 


The  presence  of  a  boundary  adjacent  to  the  grid  block 

requires  the  modification  of  equations  (3.3)  to  (3.8).  If  an 

insulator  is  in  the  ( i  -  1  ,  j  )  grid  block,  then  the  boundary 

dp 

condition  in  the  left  direction  is  *0»  This  is 
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Figure  3.1  Finite  difference  grid  used  in  MEGAERA  showing 
the  locations  at  which  potential,  temperature  and 
conductivities  are  computed. 
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approximated  by  setting  equal  to  Equation  (3.3) 

becomes 


(3.9) 


If  the  (i-1,j)  grid  block  is  part  of  an  electrode  then  0  ■  • 
is  Known.  The  coefficent  of  must  be  modified  to  take 
into  account  the  infinite  conductivity  of  the  electrode 
(which  extends  up  to  the  grid  block  boundary).  Thus,  the 
difference  equations  are 


c'j  -  0!j/C*Xi  ( h;  +  h„,  )]  (3.10) 

A',y  -  -  R,j  -  5,j  -  C'j -F, (3.11) 


(3.12) 


When  a  medium  of  different  conductivity  is  in  the 
(i-1,j)  grid  block  the  equations  (3.3)  to  (3.8)  must  be 
modified  to  deal  explicitly  with  the  conductivity 
transition.  Equations  (2.20)  and  (2.21)  give  the  boundary 
conditions  at  the  interface.  If  the  electric  field  is 
constant  between  the  center  of  the  grid  block  and  the 

t 

interface,  the  continuity  condition  for  the  electrical 
potential  (Equ.  2.21)  yields 


(3.13) 


Figure  (3.2)  shows  a  profile  of  ^  vs  x  between  the  grid 
block  centers.  Writing  equation  (2.20)  as 
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Figure  3.2  Profile  of  Electrical  Potential  S4-'  Across  an 
Conductivity  Discontinuity,  showing  the  abrupt  change  in 
normal  component  of  the  electric  field  as  required  for 
continuity  of  the  normal  component  of  current 
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(3.14) 


the  derivative 


(3.13)  may  be  eliminated. 


(3.15) 


Thus,  the  coeffi  cents  C;j  and  A  (J  in  the  finite 
difference  equation  (3.3)  become 


4ory/Ef/v  +  /v,)r*X; 


(3.16) 


(3.17) 


3.2  Differencing  the  Thermal  Equation 

The  thermal  partial  differential  equation  (2.24)  is  in 
many  respects  analogous  to  the  electrical  current  continuity 
equation  (2.10).  The  differencing  of  the  term 


is  directly  analogous  to  the  differencing  of  the  x  term  in 
equation  (2.10).  The  handling  of  the  boundary  conditions  is 
also  analogous,  with  constant  temperature  boundaries 
corresponding  to  electrodes  and  thermal  insulators 
corresponding  to  electrical  insulators. 

The  time  derivative  of  the  thermal  equation  is 
approximated  by  the  forward  difference 


(3.19) 
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where  the  superscript  indicates  the  number  of  the  timestep. 
The  heating  rate  term  is  also  calculated  and  included  in  the 
complete  difference  equation. 


3.3  Solution  of  the  Difference  Equations 

In  solving  a  set  of  partial  differential  equations  by 
the  finite  difference  technique  it  is  desirable  to  reduce 
the  set  of  equations  by  eliminating  all  but  one  unknown  (one 
unknown  being  a  single  variable  unknown  at  all  points  in  the 
domain).  This  may  be  done  before  or  after  differencing  the 
equations.  In  the  electrical  heating  simulation  presented  in 
this  thesis,  the  electrical  and  thermal  equations  are 
coupled  non- linearly  making  the  elimination  of  either  the 
potential  or  the  temperature  difficult.  For  this  reason  the 
equations  were  coupled  explicitly:  first  the  electric  field 
was  found,  then  the  new  temperatures  were  found.  This 
explicit  coupling  leads  to  convergence  problems  which  are 
discussed  in  the  section  on  program  testing. 

The  procedure  followed  by  MEGAERA  for  each  time  step  is 
as  follows: 

1.  The  electrical  conductivity  of  each  grid  block  is 
updated  from  the  temperature  of  the  grid  block  at  the 
end  of  the  last  time  step. 

2.  The  coefficents  of  the  electrical  difference  equation 
are  found  for  each  grid  block. 
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3.  This  set  of  electrical  difference  equations  is 
simultaneous.  With  one  equation  per  grid  block,  a  matrix 
of  2500  rows  and  columns  is  formed  by  MEGAERA  if  the 
full  fifty  by  fifty  grid  is  used.  This  matrix  is  solved 
by  the  Alternating  Direction  Implicit  Procedure.  This  is 
an  iterative  relaxation  method,  which  requires  an 
initial  guess  at  the  solution.  As  the  potential  does  not 
vary  greatly  from  time  step  to  time  step,  the  potential 
found  in  the  last  time  step  is  used  as  the  initial  guess 
for  the  potential  of  this  time  step  and  the  new 
potential  is  found  after  only  a  few  iterations. 

4.  With  the  potential  now  known,  the  volumetric  heating 

rate  in  each  grid  block  is  found.  MEGAERA  calculates  the 
heating  rate  directly  from  the  potential  and  cannot 
output  the  electric  field.  The  total  current  passing 
through  the  formation  between  the  electrodes  is 

calculated.  If  constant  current  or  constant  power 

heating  is  requested,  the  heating  rates  are  scaled 

* 

appropr i ately . 

5.  The  temperatures  at  the  end  of  the  time  step  are  found 
by  solving  the' thermal  equation  using  the  Alternating 
Direction  Implicit  method.  Despite  the  similarity  of  the 
names,  the  ADI  method  is  different  then  the  Alternating 
Direction  Implicit  Procedure,  the  former  being  used  for 
parabolic  equations  and  the  latter  for  elliptic 
equations . 

An  energy  balance  is  done  for  the  time  step.  The  energy 


6. 
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input,  both  for  the  time  step  and  cumulative,  is 
calculated  in  three  different  ways,  as  will  be  discussed 
in  section  3.4.1.  The  energy  calculations  are  done  for 
each  region  of  the  domain  separately,  allowing  one  to 
find  the  percentage  of  energy  input  to,  for  example,  the 
overburden  and  comparing  it  to  the  energy  input  to  the 
oi 1  sand . 

MEGAERA  was  written  to  allow  variable  grid  widths.  The 
ratio  of  maximum  width  to  minimum  width  should  not  excede 
ten  in  the  domain  of  the  problem,  with  adjacent  grid  widths 
being  within  a  factor  of  two.  The  program  is  only  two 
dimensional  and  the  maximum  number  of  grid  blocks  in  each 
the  X  and  Y  directions  is  fifty.  MEGAERA  was  written  in 
single  precision,  which  on  the  Amdahl  470  is  32  bit, 
hexidecimal  format  floating  point  arithmetic.  Using  the  full 
fifty  by  fifty  grid,  MEGAERA  uses  approximately  1.9  seconds 
of  C.P.U.  time  per  time  step  on  an  Amdahl  470 V / 7 ,  and 
requires  about  110  pages  of  virtual  memory.  Thus,  a  typical 
eighty  time  step  run  will  require  about  two  and  a  half 
minutes  of  C.P.U.  time. 


3.4  Testing  of  MEGAERA 

A  large  simulation  code  may  be  tested  in  several  ways. 
Internal  checks,  such  as  energy  balance  calculations,  give 
one  indication  of  the  accuracy  of  the  program.  Analytic 
solution  of  simple  problems  may  be  possible,  and  the  results 
of  simulating  these  problems  may  be  compared  to  the  analytic 
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solution.  Physical  modeling  of  more  complex  problems  may  be 
possible  and  numerical  simulation  results  may  be  compared  to 
physical  model  results.  Finally,  one  may  check  that 
alteration  of  the  problem  grid  and  scale  do  not 
significantly  affect  the  results  produced  by  the  simulator. 
All  of  these  methods  were  used  to  test  MEGAERA. 


3.4.1  Energy  Ba 1 ance 

The  cumulative  electrical  energy  input  at  the 
electrodes  is  found  from 


,n 


e;  r  e*'1  +  vnr*  n 


n  -r*  .  4.n 


(3.20) 


where  denotes  the  total  electrical  energy  input  up  till 

the  end  of  the  nth  time  step,  V°  and  ln  are  the  voltage 
and  current  respectively  during  the  nth  time  step,  and 
indicates  the  length  in  seconds  of  the  nth  time  step. 

The  cumulative  electrical  energy  that  is  converted  to 
heat  is  found  by  integrating  the  volumetric  heating  rate 
over  the  volume  and  time: 


(3.21 ) 


where  E#  is  the  cumulative  energy  converted  to  heat  at  the 
end  of  the  nth  time  step,  W.j  is  the  volume  of  (i,j)  grid 
block  and  Ol ..is  the  volumetric  heating  rate  in  the  grid 
block  during  the  nth  time  step. 

The  total  heat  energy  stored  in  the  temperature  rise  of 
the  formation  is  found  by  summing  the  product  of  temperature 
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rise,  volumetric  heat  capacity  and  volume  of  each  grid 
b 1 ocK : 


All  three  calculations  of  the  energy  input  to  the 
formation  should  yield  the  same  quantitative  result.  The 
computer  program  calculated  E£,  EH,  and  £Lt  for  each 
times  tep  that  was  output.  Comparing  E£  and  Ew,  the 
differences  between  both  the  cumulative  and  the  incremental 
values  of  these  quantities  were  under  one  percent  for  all 
production  runs  of  MEGAERA.  The  differences  were  greatest 
for  the  first  time  step,  usually  about  .5  to  .7  percent,  and 
decrease  as  the  run  progressed.  By  the  eightyth  time  step 
both  cumulative  and  incremental  differences  were  below  one 
quarter  percent.  The  difference  between  and  £t  was 
usually  less  than  0.1  percent.  These  figures  are  taken  from 
production  runs  (Run  No.s  10  -  20),  where  the  electrical 
conductivity  varied  from  one  region  to  the  next  by  up  to  a 
factor  of  ten. 

3.4.2  Comparison  with  Analytic  Results 

A  number  of  MEGAERA  runs  were  done  for  simple  one 
dimensional  problems  whose  solutions  could  be  found 
analytically.  Initially  the  temperature  dependence  of 
electrical  conductivity  was  set  equal  to  zero,  uncoupling 
the  electrical  equations  from  the  thermal  equation.  The 
problem  geometry  consisted  of  two  materials  of  different 
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conductivity  in  series.  The  results  showed  that  the 
potential  solver,  constant  current  or  power  scaling,  heating 
rate  calculator,  and  temperature  solver  were  working 
properly,  at  least  for  one  dimensional  problems.  In  these 
tests  the  calculated  values  agreed  with  the  analytic 
solutions  to  within  .1  percent.  The  problems  were  solved 
using  grids  of  varying  grid  width  and  differing  overall 
scale  (i.e.  maximum  dimensions). 

It  is  also  possible  to  analytically  solve  a  one 
dimensional  problem  involving  electrical  conductivity  which 
varies  with  temperature.  If  a  rectangle  of  homogeneous 
material  is  placed  between  two  parallel  plates  (and  the 
remaining  four  sides  are  insulated)  and  a  constant  voltage 
is  applied  between  the  plates,  the  electric  field  will 
remain  constant  between  the  plates.  As  the  conductivity 
rises  with  temperature  so  will  the  heating  rate.  The 
temperature  at  any  time  may  be  found  by  solving  the  integral 
equation 


(3.23) 


0 


where  d  is  the  distance  between  the  plates,  and  V  is  the 
applied  voltage  (r.m.s.  value).  This  integral  equation  may 
be  solved  by  use  of  a  Laplace  transform.  The  solution 
is22  23 


TfrMsfc  +(T,-24  <3-24) 


Runs  was  carried  out  modeling  this  problem.  When  the  maximum 
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allowed  temperature  change  per  timestep  was  three  percent, 
the  difference  between  the  computer  solution  and  the 
analytic  solution  was  one  half  of  one  percent  after  a 
temperature  rise  of  230° C.  If  the  maximum  allowed 
temperature  change  per  timestep  was  increased  to  ten 
percent,  this  difference,  after  a  similar  temperature  rise, 
was  five  percent.  The  former  case  required  roughly  three 
times  the  C.P.U.  usage  of  the  latter  case. 

The  difference  between  the  analytic  solution  and  the 
•computed  temperatures  is  due  to  the  explicit  coupling  of  the 
electrical  and  thermal  equations.  The  electrical  field  and 
heating  rates  are  held  fixed  in  a  time  step,  resulting  in  a 
conservative  approximation  to  the  actual  heating.  As  the 
conductivity  increases  with  temperature,  the  calculated 
heating  rate  will  always  be  less  then  its  actual  value.  The 
error  induced  by  this  problem  is  negligible  if  the  maximum 
allowed  temperature  change  per  timestep  is  five  percent  or 
less . 

3.4.3  Comparison  with  Physical  Models 

MEGAERA  was  used  to  simulate  a  physical  model  run  which 
was  carried  out  by  J .  Fearn  and  A.  Vogan  of  the  Applied 
Electromagnetics  group  at  the  University  of  Alberta.  The 
geometry  of  the  physical  model  was  a  scaled  down  version  of 
figure  4.2.  The  model  consisted  of  eighteen  inches  (457  mm) 
of  underlying  sand  covered  by  seven  and  three  quarters 
inches  (197  mm)  of  low  conductivity  oil  sand,  which  in  turn 
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is  covered  by  six  inches  (152  mm)  of  overburden  sand.  The 
spacing  between  the  electrodes  was  thirteen  inches  (330  mm) 
and  the  third  dimension  of  the  model,  i.e.  the  thickness  of 
the  two-dimensional  slice,  was  four  inches  (102  mm).  The 
electrodes  were  constructed  of  one  quarter  inch  (6  mm) 
copper  plate  and  were  two  inches  (51  mm)  high  and  their 
lengthes  were  equal  to  the  thickness  of  the  two-dimensional 
slice,  i.e.  four  inches  (102  mm).  The  electrodes  were 
positioned  so  the  right  electrode  was  forty  percent  in  the 
overburden  sand  and  the  left  electrode  was  forty  percent  in 
the  underlying  sand,  the  remainder  of  each  electrode  being 
in  the  oi 1  sand . 

The  underlying  and  overburden  formations  in  the  scale 
model  were  modeled  using  packed  sand  whose  moisture  content 
had  been  adjusted  so  the  conductivity  of  the  packed  sand 
(density  of  1.8  gm/cm3 )  was  7.4  X  10' 3  S/m.  The  oil  sand  was 
modeled  with  oil  sand  whose  conductivity  when  packed 
(density  of  2.0  gm/cm3)  was  2.0  X  10'3  S/m. 

The  scale  model  was  heated  by  passing  a  constant 
current  of  .183  A  between  the  electrodes  for  3605  seconds. 
The  initial  temperature  of  the  model  was  22*C  (±.8‘).  After 
the  run  was  completed  the  thermal  properties  of  the  sand  and 
oil  sand  used  in  the  model  were  measured.  The  thermal 
conductivity  of  the  underlying  and  overlying  sand  was  1.56 
W/K-m  and  its  volumetric  heat  capacity  was  1 . 37  X  106 
d/K-m3.  The  thermal  conductivity  of  the  oil  sand  was  1.61 
W/K-m  and  its  volumetric  heat  capacity  was  1.8  X  106  d/K-m3. 
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Figure  3.3  Computer  calculated  and  measured  temperatures  fo£ 
the  physical  model  run.  The  contour  plots  are  computer 
generated  and  the  measured  values  are  written  on  top. 
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The  simulation  was  carried  out  with  the  thermal  conducitity 
of  the  electrodes  initially  set  to  that  of  oil  sand,  and 
secondly  to  a  value  ten  times  greater.  This  change  caused 
only  a  slight  change  in  the  final  temperature  profile  in  the 
immediate  vicinity  of  the  electrodes. 

The  geometry,  dimensions,  sand  and  oil  sand  electrical 

and  thermal  properties,  and  the  heating  type  and  time  were 

input  to  MEGAERA.  The  temperature  dependence  of  electrical 

conductivity  was  not  measured,  but  was  assumed  to  be  0.0229 

K “ 1 ,  which  is  a  typical  value  at  low  frequencies7.  A  contour 

plot  of  the  simulation  results  is  given  in  figure  3.4  with 

the  measured  temperatures  from  the  physical  model  noted 

where  measured.  Agreement  between  the  simulation  results  and 

the  physical  model  results  is  excellent,  except  in  the 

immediate  vicinity  of  the  electrodes.  The  measured  and 

* 

computed  resistances  were  also  compared  and  the  percent 
difference  between  the  two  was  always  less  then  two  percent. 


4.  Results  of  MEGAERA  Simulation  Runs 
There  are  a  large  number  of  variables  which  affect  the 
electrical  heating  of  an  oil  sand  formation.  These  include 
the  electrical  and  thermal  properties  of  the  oil  sand  and 
surrounding  formations,  the  length  and  thickness  of  the 
electrodes,  the  spacing  between  electrodes,  the  thickness  of 
the  oil  sand,  and  the  rate  of  heating.  In  studying  the 
effects  of  varying  one  of  these  variables  the  other 
variables  must  be  held  fixed. 

The  initial  studies  done  with  MEGAERA  were  on  the 
effects  of  different  ratios  of  electrical  conductivity  of 
oil  sand  to  that  of  the  surrounding  formations,  and  on  the 
effect  of  changing  the  spacing  between  the  electrodes. 


4.1  Effects  of  the  Electrical  Conductivity  Ratio 

The  author  has  conducted  a  series  of  runs  of  MEGAERA  in 
which  the  physical  dimensions,  thermal  properties  and 
heating  rates  were  held  fixed,  but  the  electrical 
conductivities  of  the  oil  sand  and  surrounding  formations 
were  varied.  The  purpose  of  the  series  of  runs  was  to 
determine  which  of  two  electrode  configurations  produces  the 
most  favourable  temperature  distribution  for  a  given  ratio 
between  the  electrical  conductivity  of  the  oil  sand  and  that 
of  the  surrounding  formations.  The  problem  geometries, 
labeled  configuration  A  and  B,  are  shown  in  figures  4.1  and 
4.2.  In  all  runs  constant  power  heating  of  eight  kilowatts 
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Figure  4.1  Heating  Configuration  A. 
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Figure  4.2  Heating  Configuration  B.  The  thicknesses  of  these 
two  dimensional  geometries  are  set  to  one  meter  for 
resistance  and  power  calculations. 
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per  meter  length  of  the  parallel  plate  electrodes  was 
applied  for  a  period  of  one  year.  Thus,  all  the  runs  had  the 
same  heating  rate  regardless  of  the  actual  resistance 
between  the  electrodes.  The  results  of  this  series  of  runs 
are  summarized  in  Table  4.1  and  the  temperature  plots  at  the 
end  of  the  heating  period  are  given  in  Figures  4.3  to  4.9. 

In  configuration  A  (i.e.  the  electrodes  positioned  in 
the  middle  of  the  oil  sand  formation)  the  maximum 
temperature,  which  occurs  near  the  electrodes,  increases 
linearly  with  the  ratio  of  the  electrical  conductivities  of 
the  surrounding  formations  to  that  of  the  oil  sand.  (This 
ratio,  electrical  conductivity  of  overburden  and  underburden 
to  electrical  conductivity  of  oil  sand,  shall  hereafter  be 
refered  to  as  the  conductivity  ratio.)  The  temperature  in 
the  formation  mid  way  between  the  electrodes  dropped  as  the 
conductivity  ratio  increased.  Configuration  B  had  the 
opposite  trend  with  the  electrode  temperature  dropping  and 
the  mid  formation  temperature  rising  with  increased 
conductivity  ratio.  For  conductivity  ratios  greater  than  two 
configuration  B  had  a  more  desirable  temperature  profile 
than  configuration  A.  A  graph  of  temperature  vs. 
conductivity  ratio  is  given  in  figure  4.10. 

When  the  fraction  of  the  total  energy  input  that  is 
retained  in  the  oil  sand  formation  is  computed, 
configuration  A  is  superior  to  configuration  B  even  with  the 
conductivity  ratio  as  high  as  five.  With  a  conductivity 
ratio  of  one  half,  configuration  A  heating  will  retain  about 
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Figure  4.3  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  KW/m.  Run  no.  10,  configuration  B,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  10. 
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Figure  4.4  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  KW/m.  Run  no.  12,  configuration  B,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  2 . 
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Figure  4.5  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  kW/m.  Run  no.  18,  configuration  A,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  2 . 
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Figure  4.6  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  kW/m.  Run  no.  19,  configuration  A,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi  1  sand )  of  5 . 
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Figure  4.7  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  kW/m.  Run  no.  20,  configuration  B,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  5 . 
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Figure  4.8  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  kW/m.  Run  no.  21,  configuration  A,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  1 . 
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Figure  4.9  Contour  plot  of  temperature  in  Celsius  after  one 
year  of  heating  at  8  KW/m.  Run  no.  22,  configuration  A,  with 
an  electrical  conductivity  ratio  (overburden,  underburden/ 

oi 1  sand )  of  . 5 
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Figure  4.10  Temperature  vs.  the  ratio  of  the  electrical 
conductivities  of  the  surrounding  formations  to  that  of  the 

oi 1  sand . 
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able  4.1  Temperature  rise  and  thermal  efficency  as  a 
unction  of  the  relative  conductivity  of  the  oil  sand  and 

surrounding  formations. 

The  conductivity  of  the  oil  sand  is  1.0  X  10~1 2 3 4 5 6 7 8  S/m  for  all 
runs.  Constant  power  heating  of  8  KW  per  meter  was  applied 
for  one  year.  Initial  reservoir  temperature  is  15°C. 


1 . 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

Run 
No . 

Config¬ 

uration 

S/m 

x10-3 

S/m 

x10-3 

Max . 
Temp . 
‘C 

Mid  Heating 

T  emp .  i n 

°C  Oi  1 

Sand 

% 

Energy 

in 

Oi  1 
Sand 

% 

22 

A 

0.5 

0.5 

202 

101 

88 

81 

21 

A 

1.0 

1.0 

210 

85 

84 

78 

18 

A 

2.0 

2.0 

220 

67 

81 

75 

12 

B 

2.0 

2.0 

212 

65 

65 

62 

19 

A 

5.0 

5.0 

257 

50 

84 

75 

20 

B 

5.0 

5.0 

180 

66 

68 

65 

10 

B 

10. 

10.. 

177 

75 

74 

69 

1.  Runs  were  numbered  in  chronological  order. 

2.  Dimensions  and  electrode  positions  are  given  in  figures 
4 . 1  and  4.2. 

3.  Conductivity  of  overlying  formation. 

4.  Conductivity  of  underlying  formation. 

5.  Maximum  temperature  near  the  electrodes. 

6.  Maximum  temperature  mid  way  between  the  electrodes. 

7.  Electrical  energy  dissipated  in  the  oil  sand  formation 
as  a  percentage  of  the  total  energy  dissipated  in  the 
oil  sand  and  surrounding  formations. 

8.  The  energy  stored  in  the  oil  sand  formation  at  the  end 
of  one  year,  as  a  percentage  of  the  total  energy  input. 
The  difference  between  columns  7  and  8  is  due  to  thermal 
conduction . 
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eighty  percent  of  the  energy  in  the  oil  sand  after  heating 
for  one  year .  When  the  conductivity  ratio  is  increased  to 
five,  configuration  A  will  retain  three  quarters  of  the 
energy  in  the  oil  sand  formation,  while  heating  with 
configuration  B  wi 1 1  result  in  only  sixty  five  percent  of 
the  energy  input  remaining  in  the  oil  sand  at  the  end  of  the 
year.  In  both  configurations  most  of  the  energy  is  lost  by 
current  passing  through  (and  thus  heating)  the  overlying  and 
underlying  formations.  Thermal  conduction  typically  accounts 
for  the  loss  of  five  to  ten  percent  of  the  total  energy 
input,  or  ten  to  forty  percent  of  the  total  energy  lost  to 
the  surrounding  formations. 

When  heating  in  configuration  A,  the  conductivity  ratio 
determines  what  fraction  of  the  current  enters  the 
surrounding  formations.  As  the  conductivity  ratio  increases 

t 

more  of  the  current  travels  most  of  the  distance  between  the 
electrodes  in  the  overlying  or  underlying  formations,  where 
relatively  low  resistance  is  encountered.  As  a  result  the 
bulk  of  the  voltage  drop  and  the  heating  is  in  the  oil  sand 
near  the  electrodes,  where  the  current  is  travelling  from 
the  electrode  almost  directly  into  the  surrounding 
formations.  This  leads  to  a  higher  electrode  temperature  and 
a  lower  temperature  midway  between  the  electrodes.  When 
heating  in  configuration  B,  a  high  conductivity  ratio 
results  in  the  overlying  and  underlying  formations  acting 
like  large  extended  electrodes.  Most  of  the  voltage  drop  is 
across  the  oil  sand  formation  from  the  overlying  formation 
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to  the  underlying  formation.  However,  in  configuration  B 
most  of  the  current  must  travel  some  of  the  distance  between 
the  electrodes  through  the  surrounding  formations  and  this 
results  in  a  lower  fraction  of  the  total  heating  occur ing  in 
the  oil  sand  formation.  This  explains  why,  for  a 
conductivity  ratio  of  five,  configuration  B  has  a  more 
uniform  temperature  profile  in  the  oil  sand,  while 
configuration  A  has  more  of  the  energy  input  retained  in  the 
oi 1  sand . 

The  total  oil  sand  volume  in  both  configurations  A  and 
B,  per  meter  thick  slice  of  the  formation,  is  1500  m3 .  The 
chemical  energy  content  of  oil  sand  is  of  the  order  of  1010 
Joules/m3,  so  the  total  chemical  energy  in  the  heated 
formation  is  about  1013  Joules.  An  electrical  preheat  of  8 
kW/m  for  one  year  uses  about  2.5  X  1011  Joules  of  electrical 
energy,  which  is  of  the  order  of  one  percent  of  the  chemical 
energy  in  place. 


4.2  Effects  of  changing  the  distance  between  the  electrodes 
Computer  runs  were  carried  out  with  75  m  and  100  m 
spacings  between  the  electrodes,  with  the  other  dimensions 
and  electrode  positions  as  in  configuration  A  (figure  4.1). 
The  power  input  was  increased  to  12  kW/m  for  the  75  m 
spacing  and  to  16  kW/m  for  the  100  m  spacing,  but  the 
heating  time  remained  at  one  year.  The  results  of  these 
runs,  and  similar  runs  done  with  50  m  spacings,  are  given  in 
table  4.2.  The  temperature  contour  plots  for  the  end  of  the 
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year  of  heating  are  given  in  figures  4.11  to  4.17. 

Even  though  the  power  output  from  the  electrodes  is 
increased  in  proportion  to  the  distance  between  the 
electrodes,  the  maximum  temperature  near  the  electrodes 
remained  constant  for  conductivity  ratios  of  one  half  and 
one.  For  the  runs  with  different  electrical  conductivities 
for  the  overlying  and  underlying  formations,  the  maximum 
temperature  near  the  electrodes  increased  by  fifteen  percent 
when  the  distance  between  the  electrodes  was  increased  from 
50  m  to  100  m.  The  temperature  midway  between  the  electrodes 
decreased  with  an  increase  in  the  distance  between  the 
electrodes.  A  graph  of  the  temperature  midway  between  the 
electrodes  vs.  the  distance  between  the  electrodes  is  given 
in  figure  4.18.  The  fraction  of  the  total  energy  input  that 
is  retained  in  the  oil  sand  formation  at  the  end  of  the  year 
also  declines  with  an  increase  in  the  distance  between  the 
electrodes . 

As  the  distance  between  the  electrodes  increases,  more 
of  the  current  reaches  the  surrounding  formations,  and  less 
desirable  temperature  profiles  and  thermal  efficencies 
result.  A  tradeoff  must  be  made  between  the  cost  of  the 
extra  electrodes  for  a  closely  spaced  heating  configuration 
and  the  extra  electrical  energy  and  heating  time  required  to 
obtain  the  same  temperature  profile  for  a  larger  spacing. 
The  effect  of  changing  the  distance  between  the  electrodes 
is  more  pronounced  for  higher  conductivity  ratios.  The 
optimal  distance  between  the  electrodes  will  be  a  function 
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o 
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Figure  4.11  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  one  year  at  16  KW/m.  Run  no.  23, 
configuration  A,  with  a  conductivity  ratio  (overburden, 

underburden/  oil  sand)  of  .5 
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Figure  4.12  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  one  year  at  16  KW/m.  Run  no.  24, 
configuration  A,  with  a  conductivity  ratio  (overburden, 

underburden/  oi 1  sand )  of  1 . 
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Fiqure  4.13  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  one  year  at  16  KW/m.  Run  no.  2b, 
overburden  half  as  conductive  and  underburden  twice  as 

conductive  as  oil  sand. 
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Figure  4.14  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  pne  year  at  8  kW/m.  Run  no.  26, 
overburden  half  as  conductive  and  underburden  twice  as 

conductive  as  oil  sand. 
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Figure  4.15  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  one  year  at  12  kW/m.  Run  no.  32, 
overburden  half  as  conductive  and  underburden  twice  as 

conductive  as  oil  sand. 
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Figure  4.16  Contour  plot  of  temperature  in  Celsius  at  the 
pnH  of  heating  for  one  year  at  12  kw/m.  Kun  no. 
conftguratfon  A9  with  a  conductivity  ratio  (overburden, 

underburden/  oil  sand)  of  .b 
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Figure  4.17  Contour  plot  of  temperature  in  Celsius  at  the 
end  of  heating  for  one  year  at  12  KW/m.  Run  no.  35, 
configuration  A,,  with  a  conductivity  ratio  (overburden, 

underburden/  oil  sand)  of  1. 


TEMPERATURE  MID  WAY  BETWEEN  THE  ELECTRODES  («C) 
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*  Overlying  formation  is  half  as  conductive  as  the  oil  sand 
and  the  underlying  formation  Is  twice  as  conductive 
as  the  oil  sand. 


Figure  4.18  Temperature  midway  between  the  electrodes  vs. 
the  distance  between  the  electrodes.  Electrode  positions  are 
as  in  configuration  A,  except  for  the  mixed  conductivity 
ratio  where  the  electrodes  were  moved  up  three  meters. 
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Table  4.2  Temperature  rise  and  thermal  efficency  as  a 
function  of  distance  between  the  electrodes. 

The  conductivity  of  the  oil  sand  is  1.0  X  10’1 2 3 4 5 6 7 8  S/m  for  all 

power  heating  of  8  kW  per  meter  in  the  runs 
with  a  50  m  spacing,  12  kW  per  meter  in  the  runs  with  a  75  m 
spacing,  and  16  kW  per  meter  in  the  runs  with  a  100  m 
spacing  was  applied  for  one  year.  Formation  thicknesses  and 
electrode  sizes  and  positions  are  as  in  configuration  A 


1 . 

2. 

3. 

4. 

Run 

No. 

Spacing 

m 

&QB 

S/m 

XlO-3 

OTa 

S/m 

xlO-3 

22 

50 

0.5 

0.5 

34 

75 

0.5 

0.5 

23 

100 

0.5 

0.5 

21 

50 

1.0 

1.0 

35 

75 

1.0 

1.0 

24 

100 

1.0 

1.0 

26 

50 

0.5 

2.0 

32 

75 

0.5 

2.0 

25 

100 

0.5 

2.0 

•  ■  /  i 


5. 

6. 

7. 

8. 

Max . 

Mid 

Heat i ng 

Energy 

Temp . 

Temp . 

in 

i  n 

°C 

’C 

Oi  1 

Oi  1 

Sand 

Sand 

% 

% 

202 

101 

88 

81 

196 

82 

79 

73 

195 

72 

72 

67 

210 

85 

84 

78 

210 

60 

71 

66 

210 

48 

62 

58 

230 

90 

85 

79 

258 

74 

74 

68 

264 

46 

65 

60 

1.  Runs  were  numbered  in  chronological  order. 

2.  Distance  between  the  electrodes. 

3.  Conductivity  of  overlying  formation. 

4.  Conductivity  of  underlying  formation. 

5.  Maximum  temperature  near  the  electrodes. 

6.  Maximum  temperature  mid  way  between  the  electrodes. 

7.  Electrical  energy  dissipated  in  the  oil  sand  formation 
as  a  percentage  of  the  total  energy  dissipated  in  the 
oil  sand  and  surrounding  formations. 

8.  The  energy  stored  in  the  oil  sand  formation  at  the  end 
of  one  year,  as  a  percentage  of  the  total  energy  input. 
The  difference  between  columns  7  and  8  is  due  to  thermal 
conduction . 
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of  many  parameters,  with  two  important  factors  being  the 
cost  of  the  wells  and  electrodes  and  the  relative 
conductivities  of  the  surrounding  formations  to  the  oil 


sand . 


.  ? 


5.  Conclusions 

The  results  of  the  full  scale  simulation  runs  done  using  the 
computer  program  MEGAERA  are  encouraging.  These  results 
indicate  that,  by  using  less  than  five  percent  of  the 
chemical  energy  in  place  in  a  typical  thickness  Athabasca 
oil  sand  formation,  it  is  possible  to  raise  the  temperature 
midway  between  the  electrodes  to  80° C,  with  higher 
temperatures  near  the  electrodes.  The  viscosity  of  the 
bitumen  (Athabasca  oil  sand)  would  be  about  1000  cp  midway 
between  the  electrodes,  and  considerably  less  near  the 
electrodes.  This  is  low  enough  to  allow  an  efficent  steam 
drive  to  be  conducted. 

The  comparisons  done  between  the  analytic  solution  and 
the  computer  simulation  indicate  that  the  finite  difference 
scheme  used  was  correct.  The  convergence  error  due  to  the 
explicit  coupling  of  the  electrical  and  thermal  equations 
remains,  but,  if  the  maximum  temperature  change  per  times tep 
is  kept  less  than  five  percent,  this  error  will  not  be 
signi f icant . 

The  best  test  of  both  the  mathematical  model  and  the 
finite  difference  scheme  that  were  used  is  conducted  by 
comparing  the  simulator  results  with  a  physical  model  run. 
When  this  was  done  the  simulator  results  were  in  agreement 
with  the  measured  temperatures  from  the  physical  model  run, 
except  in  the  immediate  vicinity  of  the  electrodes.  Thus,  at 
least  for  temperatures  up  to  90* C  (above  which  accurate 
electrical  conductivity  data  is  not  available),  the 
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mathematical  model  used  will  accurately  predict  temperatures 
and  resistance  for  electrical  heating  of  solid  materials. 

There  are  several  items  which  should  be  Kept  in  mind 
when  interpreting  the  results  of  the  full  scale  simulator 
runs.  First,  the  electrical  and  thermal  properties  of  oil 
sand  have  not  been  reported  in  the  liturature  for 
temperatures  above  90*C.  In  the  simulator  the  thermal 
properties  were  assumed  independent  of  temperature,  and  the 
electrical  conductivity  was  assumed  to  increase  linearly 
with  temperature,  up  to  250°C.  Second,  the  conductivity  of 
the  oil  sand  near  the  electrodes  in  a  field  test  would  be 
altered  by  the  injection  of  brine  near  the  electrode.  These 
two  factors  may  result  in  significant  error  in  the  predicted 
temperatures  near  the  electrodes.  Finally,  the  conductivity 
of  the  oil  sand  was  made  uniform  throughout  the  formation 
for  the  simulation  runs.  This  was  by  choice  and  not  due  to 
any  limitation  of  the  simulator.  In  an  actual  oil  sand 
formation  there  is  probably  variations  in  electrical 
conductivity  with  height  in  the  formation  and  due  to  the 
presence  of  shale  breaks  and  clay  deposits. 

The  investigation  to  date  indicates  that  the  electrical 
preheat  method  may  be  developed  into  an  economical  method  of 
preparing  an  oil  sand  formation  for  a  steam  flood  or  fire 
flood.  Further  measurement  of  the  electrical  and  thermal 
properties  of  oil  sand  need  to  be  performed  at  formation 
pressures  and  temperatures  in  the  range  of  100*C  to  300* C. 
The  computer  simulator  should  be  extended  to  three 
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dimensions  so  that  simulation  of  five  spot  patterns  are 
possible  in  addition  to  the  parallel  plate  electrode  (or 
line  drive)  geometry  which  was  studied  in  this  thesis. 
Finally,  due  to  variations  in  the  conductivities  of  the  oil 
sand  and  the  surrounding  formations  no  one  electrode 
configuration  is  optimal  and  each  recovery  location  will 
have  to  be  investigated  individually  to  determine  the 
optimal  electrode  size,  positioning  and  spacing. 
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Appendix  I 


Listing  of  the  Program  MEGAERA 
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********************* 
*  * 

*  MEGAERA  * 

*  * 
********************* 


SIMULATION  OF  ELECTRICAL  CONDUCTION  HEATING 


OF  OILSAND  AND  OTHER  MATERIALS 


ALLAN  HI  EBERT ,  AUGUST  1980 

WRITTEN  FOLLOWING  THE  OLYMPUS 
PROGRAMING  SYSTEM  DEVELOPED  AT 
CULHAM  LABORITORI ES ,  ENGLAND. 


INDEX  OF  SUBPROGRAMS 


MAIN  CONTROL.  CLASS  0 


MAIN 

FORTRAN  MAIN  PROGRAM 

0. 

0 

BASIC 

INITIALIZE  BASIC  CONTROL  DATA 

0. 

1 

MODIFY 

MODIFY  BASIC  DATA  IF  REQUIRED 

0. 

2 

COTROL 

CONTROL  THE  RUN 

0. 

3 

USER 

PRINT  USER,  TIME  AND  DATE  TO  DIARY 

0. 

4 

PROLOGUE . 

CLASS  1 

LABRUN 

LABEL  THE  RUN 

1.. 

1 

CLEAR 

CLEAR  VARIABLES  AND  ARRAYS 

l‘. 

2 

PRESET 

SET  DEFAULT  VALUES 

1  . 

3 

DATA 

DEFINE  DATA  SPECIFIC  TO  RUN 

1  . 

4 

AUXVAL 

SET  AUXILIARY  VALUES 

1 . 

5 

INITAL 

DEFINE  PHYSICAL  INITIAL  CONDITIONS 

1  . 

6 

RESUME 

RESUME  FROM  PREVIOUS  RECORD 

1  . 

7 

START 

START  OR  RESTART  THE  RUN 

1  . 

8 

EPARAM 

CALCULATE  ITERATION  PARAMETERS  FOR  A 

.D. I .P.  1 . 

9 

TCOEFF 

CALCULATE  CONSTANTS  FOR  THE  THERMAL 

EQUATION  1. 

10 

CALCULATION. 

CLASS  2 

STEPON 

STEP  ON  THE  CALCULATION 

2. 

1 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  10 

1  1  1 

1  12 

1  13 

1  14 

1  15 

1  16 

117 

118 

1  19 

120 
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ELEPOT 

ELECOF 

THOMAS 

MAXDIF 

ECOND 

OCALC 

TCALC 

ENGBAL 

ELECUR 


OUTPUT ( 1 ) 
OUTGRD( 1 ) 
ALINTP( 7 ) 
OUT I NT ( 2 ) 
OUTR 
OUT  I 
OUTH 

OUTTAP( 1 ) 


TESEND 

ENDRUN 


CALCULATE  THE  ELECTRICAL  POTENTIAL 
CALCULATE  COEFFICENTS  OF  DIFFERENCE  EOUATION 
SOLVE  TRIDIAGONAL  SYSTEM  BY  THOMAS  ALGORITHM 
FIND  THE  MAXIMUM  DIFFERENCE  BETWEEN  TEM3 ,  TEM1 
CALCULATE  THE  ELECTRICAL  CONDUCTIVITY 
CALCULATE  THE  HEATING  RATES  FOR  THE  GRID 
SOLVE  THE  HEAT  DIFFUSION  EOUATION 
CALCULATE  THE  ENERGY  BALANCE  FOR  THE  TIMESTEP 
CALCULATE  THE  CURRENT  THROUGH  A  SURFACE 


OUTPUT.  CLASS  3 

CONTROL  THE  OUTPUT  3  ^ 

OUTPUT  ONE  OF  THE  GLOBAL  VARIABLE  ARRAYS  3  2 

DO  A  LINEAR  INTERPOLATION  3.3 

INTERPOLATE  AND  PRINT  A  GOBAL  VARIABLE  3^4 

OUTPUT  A  LABEL  AND  A  REAL  VARIABLE  3.5 

OUTPUT  A  LABEL  AND  AN  INTEGER  VARIABLE  3.6 

OUTPUT  A  LABEL  AND  A  HOLLERITH  VARIABLE  3.7 

OUTPUT  RUN  INFORMATION  TO  MAGNETIC  TAPE  3.8 

EPILOGUE.  CLASS  4 


TEST  FOR  COMPLETION  OF  RUN 
TERMINATE  THE  RUN 


2.2 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 
2.9 
2.  10 


DIAGNOSTICS.  CLASS  5 

REPORT ( 3  )  CONTROL  THE  DIAGNOSTICS  5.1 
CLIST ( 2 )  PRINT  COMMON  VARIABLES  5.2 
ARRAYS ( 2  )  PRINT  COMMON  ARRAYS  5.3 


UTILITIES.  CLASS  U 

MESAGE ( 1  )  PRINT  48-CHARACTER  MESSAGE  ON  OUTPUT  CHANNEL  U.1 
PAGE  FETCH  NEW  PAGE  ON  OUTPUT  CHANNEL  U.2 

BLINES(I)  INSERT  BLANK  LINES  ON  OUTPUT  CHANNEL  U.3 

RVAR  ( 2 )  PRINT  NAME  AND  VALUE  OF  REA*L  VARIABLE  U.4 

I  VAR ( 2 )  PRINT  NAME  AND  VALUE  OF  INTEGER  VARIABLE  U.5 

HVAR ( 2 )  PRINT  NAME  AND  VALUE  OF  HOLLERITH  VARIABLE  U.6 

LVAR ( 2 )  PRINT  NAME  AND  VALUE  OF  LOGICAL  VARIABLE  U.7 

RARRA Y ( 3 )  PRINT  NAME  AND  VALUES  OF  REAL  ARRAY  U.8 

I  ARRAY ( 3  )  PRINT  NAME  AND  VALUES  OF  INTEGER  ARRAY  U.9 

HARRA Y ( 3 )  PRINT  NAME  AND  VALUES  OF  HOLLERITH  ARRAY  U.10 

REPTHDO)  PRINT  HEADING  FOR  DIAGNOSTIC  REPORT  U.11 

RUNTIM  UPDATE  CPU  TIME  AND  PRINT  IT  U.12 

DAYTIM  PRINT  DATE  AND  TIME  U.13 

RESETR( 3 )  RESET  REAL  ARRAY  TO  SPECIFIED  VALUE  U.14 

RESETI ( 3 )  RESET  INTEGER  ARRAY  TO  SPECIFIED  VALUE  U.15 

RESETH( 3 )  RESET  HOLLERITH  ARRAY  TO  SPECIFIED  VALUE  U.16 

JOBTIM  FETCH  ALLOCATED  JOBTIME  U.17 

LARRAY ( 3 )  PRINT  NAME  AND  VALUES  OF  LOGICAL  ARRAY  U.18 

RESETL( 3 )  RESET  LOGICAL  ARRAY  TO  SPECIFIED  VALUE  U.19 

RARAY2  PRINT  DOUBLY-SUBSCRIPTED  ARRAY  U . 20 

SCALER ( 3 )  SCALE  A  REAL  ARRAY  BY  A  REAL  VALUE  U.21 

SCALEI ( 3 )  SCALE  AN  INTEGER  ARRAY  BY  AN  INTEGER  VALUE  U . 22 

COPYR ( 5 )  COPY  ONE  REAL  MATRIX  INTO  ANOTHER  U . 23 

COPY  I ( 5  )  COPY  ONE  INTEGER  MATRIX  INTO  ANOTHER  U . 24 

S I GNR ( 2  )  CHANGE  THE  SIGN  OF  A  REAL  MATRIX  U . 25 

SIGNI ( 2 )  CHANGE  THE  SIGN  OF  AN  INTEGER  MATRIX  U . 26 


121 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

*157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
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DUMC0M(3)  DUMP 


SELECTED  COMMON  BLOCKS 


U .  27 


CL  Cl . 1 .  BASIC  SYSTEM  PARAMETERS 

C  VERSION  2B  14.8.73  KVR/MHH  CULHAM 

C  COMMON/COMBAS/ 

c 

c 


c 

ALTIME 

TIME  ALLOCATED  TO  JOB,  (MINS) 

R 

1  .  1 

c 

CPTIME 

CPU  TIME  USED  SO  FAR  ON  THIS  JOB,  (MINS) 

R 

1  .  1 

c 

LABEL  1 ( 12) 

LABEL 

DESCRIBING  THE  RUN 

IA 

1  .  1 

c 

LABEL2 ( 12) 

LABEL 

DESCRIBING  THE  RUN 

IA 

1  .  1 

c 

LABEL3( 12) 

LABEL 

DESCRIBING  THE  RUN 

IA 

1  .  1 

c 

LABEL4 ( 12) 

LABEL 

DESCRIBING  THE  RUN 

IA 

1  .  1 

c 

LABEL5( 12) 

LABEL 

AVAILABLE  TO  PROGRAMMER 

IA 

1  .  1 

c 

LABEL6( 12) 

LABEL 

AVAILABLE  TO  PROGRAMMER 

IA 

1  .  1 

c 

LABEL7 ( 12) 

•LABEL 

RESERVED  FOR  SYSTEM  USE 

IA 

1  .  1 

c 

LABEL8( 12) 

LABEL 

RESERVED  FOR  SYSTEM  USE 

IA 

1  .  1 

c 

NDIARY 

CHANNEL  FOR  DIARY 

I 

1  .  1 

c 

NIN 

CURRENT  INPUT  CHANNEL 

I 

1  .  1 

c 

NLEDGE 

CHANNEL  FOR  RESTART  RECORDS 

I 

1  .  1 

c 

NLEND 

.TRUE . 

IF  RUN  TO  BE  TERMINATED 

L 

1  .  1 

c 

NLRES 

.TRUE . 

IF  RUN  TO  BE  RESTARTED 

L 

1  .  1 

c 

NONLIN 

CHANNEL  FOR  ONLINE  INPUT-OUTPUT 

I 

1  .  1 

c 

NOUT 

CURRENT  OUTPUT  CHANNEL 

I 

1  .  1 

c 

NPRINT ‘ 

CHANNEL  FOR  PRINTED  OUTPUT 

I 

1  .  1 

c 

NPUNCH 

CHANNEL  FOR  CARD  OUTPUT  (OR  EQUIVALENT) 

I 

1  .  1 

c 

NREAD 

CHANNEL  FOR  CARD  INPUT 

I 

1  .  1 

c 

NREC 

CURRENT  RECORD  NUMBER 

I 

1  .  1 

c 

NRESUM 

RESUME 

FROM  RECORD  ON  THIS  CHANNEL 

I 

1  .  1 

c 

NRUN 

MAXIMUM  NUMBER  OF  STEPS 

I 

1  .  1 

c 

NSTEP 

CURRENT  STEP  NUMBER 

I 

1  .  1 

c 

STIME 

START 

TIME,  (MINS) 

R 

1  .  1 

c 

c 

c - 

c 

C 

CL  Cl. 9.  DEVELOPMENT  AND  DIAGNOSTIC  PARAMETERS 

C  VERSION  IB  14.8.73  KVR/MHH  CULHAM 

C  COMMON/ COMDDP/ 

c 

c 

c 


C 

MAXDUM 

MAXIMUM  DIMENSION  OF  DUMP  ARRAYS 

I 

1  .9 

C 

MXDUMP 

ACTUAL  DIMENSION  OF  DUMP  ARRAYS 

I 

1 .9 

C 

NADUMP(M) 

CODES  FOR  ARRAY  DUMPS 

IA 

1 .9 

C 

NCLASS 

MOST  RECENT  CLASS  REPORTED 

I 

1 .9 

C 

NLCHED 

.TRUE.  IF  CLASS  0  REPORT  HEADS  REQUIRED 

L 

1 .9 

C 

NLHEAD ( 9 ) 

.TRUE.  IF  CLASSES  1-9  REPORT  HEADS  REQUIRED 

LA 

1 .9 

C 

NL0MT1 (50) 

CLASS  1  SUBPROGRAM  SELECTOR 

LA 

1 .9 

C 

NL0MT2 ( 50 ) 

CLASS  2  SUBPROGRAM  SELECTOR 

LA 

1 .9 

C 

NL0MT3 ( 50 ) 

CLASS  3  SUBPROGRAM  SELECTOR 

LA 

1 .9 

C 

NLREPT 

.TRUE.  IF  ANY  REPORT  REQUIRED 

L 

1 .9 

C 

NPDUMP ( M ) 

CODES  FOR  DUMPING  POINTS 

IA 

1 .9 

c 

NPOINT 

MOST  RECENT  POINT  REPORTED 

I 

1 .9 

c 

NSUB 

MOST  RECENT  SUBPROGRAM  REPORTED 

I 

1 .9 

c 

NVDUMP(M) 

CODES  FOR  DUMPING  VARIABLES 

IA 

1 .9 

c 

c 


N 


' 


181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 

218 

219 

220 

221 

222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 
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c 

c 

c 


1 .2 


/C0MGL0/  -  GLOBAL  VARIABLES 


c 

DELTAX 

WIDTHS  OF  GRID  BLOCKS 

R 

1 .2 

c 

DELTAY 

HEIGHT  OF  GRID  BLOCKS 

R 

1  .2 

c 

DTIME 

LENGTH  OF  TIMESTEP  (  DELTA  TIME) 

R 

1.2 

c 

ELECON 

ELECTRICAL  CONDUCTIVITY 

R 

1 .2 

c 

ELEGEO 

ELECTRICAL  GEOMETRY 

R 

1 .2 

c 

NEGATIVE  -  DOMAIN  POINT 

c 

ZERO  -  NOT  IN  DOMAIN 

c 

POSITIVE  -  ELECTRODE 

c 

OTHERM 

HEAT  GENERATED  BY  ELECTRIC  FIELD 

R 

1 .2 

c 

POTENT 

ELECTRICAL  POTENTIAL 

R 

1 .2 

c 

TEMP 

TEMPERATURE 

R 

1 .2 

c 

THICK 

THICKNESS  OF  TWO  DIMENSIONAL  SLICE 

R 

1 .2 

c 

THMGEO 

THERMAL  GEOMETRY 

R 

1 .2 

c 

TIME 

CURRENT  VALUE  OF  THE  TIME 

R 

1.2 

c 

XCOORD 

X  COORDINATES  OF  GRIDPOINTS 

R 

1 .2 

« 

YCOORD 

Y  COORDINATES  OF  GRIDPOINTS 

R 

1 .2 

c 

ORIGIN  AT  LOWER  LEFT  CORNER  OF  (1,1) 

BLOCK 

c 

XMIN 

MINIMUM  X  FOR  INTERPOLATED  PRINTOUT 

R 

1 .2 

c 

XMAX 

MAXIMUM  X  FOR  INTERPOLATED  PRINTOUT 

R 

1 .2 

c 

YMIN 

MINIMUM  Y  FOR  INTERPOLATED  PRINTOUT 

R 

1 .2 

c 

YMAX 

MAXIMUM  Y  FOR  INTERPOLATED  PRINTOUT 

R 

1 .2 

C 

C- 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

c 


1 .3 


/COMELE/  -  ELECTRICAL  COEFFICENTS 


ECXM 

ECXP 

ECYM 

ECYP 

EXMXP 

EYMYP 

ERHS 

ELEVOL 

ELEALP 

ELEBET 

ELEPAR 

EPSELE 


COEFFICENT 

COEFFICENT 


SUM  OF  THE 
SUM  OF  THE 
RIGHT  HAND 
VOLTAGES  OF 
COEFFICENTS 


OF  FINITE  DIFFERENCE  EOU . 

IN  POSITIVE  X  DIRECTION 
"  NEGATIVE  Y  " 

"  POSITIVE  Y  DIRECTION 
TWO  X  DIRECTION  COEFFICENTS 
TWO  Y  DIRECTION  COEFFICENTS 
SIDE  OF  DIFFERENCE  EOU. 
DIFFERENT  ELECTRODES 
FOR  CONDUCTIVITY  CALC. 


ITERATION  PARAMETERS  FOR  A.D.I.P. 
CONVERGENCE  CRITERIA  FOR  POTENTIAL 


R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 

R 


1 .3 


1 .4 


/COMTHM/  -  THERMAL  CONDUCTIVITY  AND  CAPACITY 


THMCON 

THERMAL  CONDUCTIVITIES 

R 

1 . 4 

THMCAP 

VOLUMETRIC  HEAT  CAPACITY 

R 

1 . 4 

THMTEM 

VECTOR  OF  CONSTANT  TEMPERATURES 

R 

1  .  4 

TCYP 

COEFFICENT  OF  THE  THERMAL 

EOU 

.  ( 

Y 

POS. 

) 

R 

1 . 4 

TCXP 

ft  MW  M 

N 

( 

X 

POS  . 

) 

R 

1 .4 

TCXM 

M  n  M  M 

H 

( 

X 

NEG. 

) 

R 

1 .4 

TCYM 

M  M  M  M 

W 

( 

Y 

NEG. 

) 

R 

1  .  4 

TCX 

SUM  OF  X  COEFFICENTS 

R 

1 . 4 

TCY 

SUM  OF  Y  COEFFICENTS 

R 

1 . 4 

TCRHS 

TERMS  ON  CONSTANT  TEMPERATURE 

BOUNDARY 

R 

1  .  4 

1.5  /COMTEM/  -  TEMPERARY  ARRAYS 


241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 

281 

282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 
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c 

TEM1 

TWO 

DIMENSION 

ARRAY 

FOR 

c 

TEM2 

c 

TEM3 

c 

T  1 

ONE 

DIMENSION 

ARRAY 

FOR 

C 

c 

c 

c 

c 

c- 

c 

c 

c 


T2 
T3 
T  4 
T5 


R 

R 

R 

R 

R 

R 

R 

R 


1.5 
1.5 
1 .5 
1  .5 
1 .5 
1  .5 
1 .5 
1.5 


C 

C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


1.6  /COMDIM/  -  INTEGER  COMMON  BLOCK 


C 

NX 

ACTUAL 

NUMBER  OF  GRID  BLOCKS  IN  X  DIR. 

I 

1 .6 

C 

NY 

ACTUAL 

NUMBER  OF  GRID  BLOCKS  IN  Y  DIR. 

I 

1 .6 

C 

NDX 

ACTUAL 

DIMENSION  OF  ARRAYS 

I 

1 . 6 

C 

NDY 

ACTUAL 

DIMENSION  OF  ARRAYS 

I 

1 .6 

C 

NEOU 

NUMBER 

OF  EQUATION  FOR  THOMAS  ALGOR. 

I 

1 .6 

C 

NALPHA 

NUMBER 

OF  ITERATION  PARAMETERS  USED 

I 

1 .6 

C 

MITER 

CURRENT  ITERATION  COUNTER  FOR  POTENTIAL 

I 

1  .6 

C 

MM  AX 

MAXIMUM  ALLOWED  VALUE  OF  MITER 

I 

1 .6 

C 

NREG 

NUMBER 

OF  DIFFERENT  REGIONS  IN  PROBLEM 

I 

1 .6 

c 

NGEO 

0  IF  IN  VERTICAL  MODE,  1  IN  HORIZONTAL 

I 

1 .6 

c 

NERROR 

ERROR  RETURN  CODE  FROM  TlMESTEP  ROUTINES 

I 

1 .6 

c 

NTYPE 

INDICATES  TYPE  OF  HEATING 

I 

1 .6 

c 

1- 

CONSTANT  VOLTAGE 

c 

2- 

CONSTANT  CURRENT 

c 

3- 

CONSTANT  POWER 

c 

4- 

NO  HEATING 

c 

NSPO 

NUMBER 

OF  TIME  STEPS  PER  OUTPUT 

I 

1 .6 

c 

NSTO 

NUMBER 

OF  TIMESTEPS  BETWEEN  TAPE  STORES 

I 

1 .6 

c 

NPX 

NUMBER 

OF  X  INTERVALS  FOR  INTERPOLATION 

I 

1 .6 

c 

NPY 

NUMBER 

OF  Y  INTERVALS  FOR  INTERPOLATION 

I 

1 .6 

c 

NIU 

INDICATES  IF  SURFACE  OF  CURRENT  INTEGERAL 

c 

IS  PERPENDICULAR  TO  X  OR  Y 

I 

1 .6 

c 

NCI 

INDICE 

OF  GRIDLINE  WHICH  INTEGRAL  IS  ALONG 

I 

1 .6 

1  .7 


/COMCON/  -  HEATING  CONTROL  VARIABLES 


HTIME 

TOTAL  TIME  OF  HEATING 

R 

1 .7 

CTIME 

TOTOL  TIME  OF  COOLING 

R 

1 . 7 

CCUR 

CONSTANT  VALUE  OF  CURRENT 

R 

1 .7 

CPOW 

CONSTANT  VALUE  OF  CURRENT 

R 

1 . 7 

DTEMP 

MAXIMUM  RELATIVE  CHANGE  IN  TEMPERATURE 

R 

1  .  7 

DELT 

ACTUAL  CHANGE  IN  TEMPERATURE  FOR  TlMESTEP 

R 

1 . 7 

CUR 

CURRENT  CALCULATED  FORM  CURRENT  INTEGRAL 

R 

1 .7 

VOLTS 

ACTUAL  VOLTAGE  ACROSS  ELECTRODES 

R 

1  .  7 

CVOL 

VOLTAGE  USED  IN  PROBLEM  DEFINITION 

R 

1 . 7 

RESIST 

resistance  BETWEEN  ELECTRODES 

R 

1  .  7 

POWER 

ELECTRICAL  POWER  INPUT  DURING  TlMESTEP 

R 

1  .  7 

TINIT 

INITIAL  TEMPERATURE  OF  THE  FORMATION 

R 

1  .7 

DEENG 

ELECTRICAL  ENERGY  INPUT  FOR  THE  TlMESTEP 

R 

1 . 7 

TEENG 

TOTAL  ELECTRICAL  ENERGY  INPUT  FOR  THE  RUN 

R 

1 . 7 

DQENG 

HEAT  PRODUCED  FOR  THE  TlMESTEP 

R 

1 . 7 

TQENG 

TOTAL  HEAT  PRODUCED  IN  THE  RUN 

R 

1 . 7 

QENG 

VECTOR  OF  HEATING  IN  EACH  REGION 

R 

1 . 7 

TTENG 

ENERGY  IN  THE  TEMPERATURE  CHANGE 

R 

1 . 7 

V  i  * 


301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 
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C  TENG  VECTOR  OF  ENERGY  OF  TEMP.  CHANGE  BY  REGION  R  1.7 

C  SCALE  SCALE  FACTOR  FOR  CONSTANT  CURRENT  AND  POWER  R  1.7 

C 

C - 

c 

c 

C  0.0  FORTRAN  MAIN  PROGRAM 
C 

c 

C/  INSERT  COMBAS 
C 

C  TIME  ALLOCATED  TO  JOB 

CALL  JOBTIM(ALTIME) 

C 

C  SET  UP  THE  BASIC  CONTROL  DATA 

CALL  BASIC 
C 

C  PRINT  DATE  AND  TIME 

CALL  PAGE 
CALL  DAYTIM 
C 

C  CONTROL  THE  RUN 

CALL  COTROL 
C 

STOP 

END 

C 

C 

- - 

c 

SUBROUTINE  BASIC 
C 

C  0.1  INITIALIZE  BASIC  DATA 
C 

c 

C/  INSERT  COMBAS 
C/  INSERT  COMDDP 

DATA  IBLANK/4H  / 

CL  1 .  GENERAL  OLYMPUS  DATA 

C 

CL  i.i  BASIC  SYSTEM  PARAMETERS 

C  CPU  -  TIME  USED  SO  FAR 

CPTIME  =0 . 0 
C 

c  CLEAR  ALL  8  LABEL  ARRAYS 

IL  =  8  *  12 

CALL  RE SETH ( LABEL  1 , I L , I BLANK ) 

C 

C  INPUT-OUTPUT  CHANNELS 

NLEDGE  =  30 
NONLIN  =  1 
NPUNCH=7 
NPRINT=6 
NREAD=5 
NDI ARY  =NPR INT 
NIN=NREAD 
NOUT  =  NPRINT 
C 

C  TIMESTEP  CONTROL 

NRUN  *  1 


"N 
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367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 

408 

409 

410 

411 

412 
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415 

416 

417 
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419 
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NSTEP=0 

C 

c  RESTART  CONTROL 

NREC  =  1 
NRESUM  =  NLEDGE 
C 

C  LOGICAL  SWITCHES 

NLEND= . FALSE . 

NLRES= . FALSE . 

C 

CL  1-9  DIAGNOSTIC  AND  DEVELOPMENT  PARAMETERS 

MAXDUM  =  20 

C  MAXIMUM  DIMENSIONS  OF  DUMP  ARRAYS 

MXDUMP  =  10 

•  C  RESET  DUMP  ARRAYS 

CALL  RESETI (NADUMP , MAXDUM, 0) 

CALL  RESETI (NPDUMP . MAXDUM, 0) 

CALL  RESETI (NVDUMP , MAXDUM, 0) 

C  TRACER  VARIABLES 

NCLASS  =  0 
NSUB  =  1 
NPOINT  =  1 

C  LOGICAL  SWITCHES 

NLCHED  =  .FALSE. 

NLREPT  =  .FALSE. 

C  REPORT  HEADS  FOR  CLASSES  1-9 

CALL  RESETL(NLHEAD,9, .FALSE. ) 

C  RESET  CLASS  1,2,3  SUBPROGRAM  SELECTOR  ARRAY 

CALL  RESETL(NL0MT1 ,50, .FALSE . ) 

CALL  RESETL(NL0MT2, 50,  .FALSE.  ) 

CALL  RESETL(NL0MT3, 50,  .FALSE.  ) 

C 

C  USER  INTERFACE 

CALL  MODIFY 
C 

RETURN 

END 

C 

c - 

c 

SUBROUTINE  MODIFY 
C 

C  0.2  MODIFY  BASIC  DATA  IF  REQUIRED 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMDDP 
C 

NAMELIST / INMOD/NONL IN , NRUN , NREC , NRESUM , NLRES 
C 

C  PRINT  USER,  TIME  AND  DATE  TO  DIARY  FILE 

CALL  USER 
C 

C  SET  DEFAULT  VALUE  FOR  NRESUM 

NRESUM=2 
C 

C  INPUT  NAMELIST  TO  MODIFY  BASIC  DATA 

READ ( NREAD , INMOD) 

I F ( .NOT. NLRES)  NREC=1 
RETURN 
END 
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C 

- - 

c 

SUBROUTINE  COTROL 


c 

C  0.3 

CONTROL  THE  RUN 

C 

C 

C 

VERSION  2B 

17/12/73  KVR/MHH  CULHAM 

C/  INSERT  COMBAS 

C/  INSERT  COMDDP 

CL 

1  . 

PROLOGUE 

C 

IF(NLRES)  GO  TO 

170 

C 

C 

A  . 

NEW  RUN 

c 

CL 

1  .  1 

LABEL  THE  RUN 

1  10 

CALL  LABRUN 

C 

CL 

1  .  2 

CLEAR  VARIABLES  AND  ARRAYS 

120 

CALL  CLEAR 

C 

CL 

1 .3 

SET  DEFAULT  VALUES 

130 

CALL  PRESET 

C 

CL 

1 .4 

DEFINE  DATA  SPECIFIC  TO  RUN 

140 

CALL  DATA 

C 

CL 

1  .5 

SET  AUXILIARY  VALUES 

150 

CALL  AUXVAL 

C 

CL 

1 .6 

DEFINE  PHYSICAL  INITIAL  CONDITIONS 

160 

CALL  INITAL 

C 

C 

GO  TO  180 

B  . 

RESUME  A  PREVIOUS  RUN 

C 

CL 

1 . 7 

PICK  UP  RECORD,  MODIFY  REQUIRED  PARAMETERS 

170 

CONTINUE 

C 

LABEL  THE  CONTINUATION  RUN 

C 

CALL  LABRUN 

CLEAR  VARIABLES  AND 

ARRAYS 

C 

CALL  CLEAR 

PICK  UP  RECORD  AND 

PRINT  DETAILS 

C 

CALL  RESUME 

READ  ANY  NEW  DATA  NEEDED 

C 

CALL  DATA 

MODIFY  AUXILIARY  VARIABLES  AS  REQUIRED 

C 

c 

CALL  AUXVAL 

C  . 

PRELIMINARY  OPERATIONS 

c 

CL 

1  .8 

START  OR  RESTART  THE  RUN 

180 

CALL  START 

C 

INITIAL  OUTPUT 

C 

C 

CL 

CALL  OUTPUT ( 1 ) 

2  . 

CALCULATION 

' 
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c 

CL 

210 
C 

CL 
C 

CL 

310 
C 

CL 
C 

CL 

410 
C 

C  FINAL  OUTPUT 

CALL  OUTPUT ( 3 ) 

C 

CL  4.2  TERMINATE  THE  RUN 

420  CALL  ENDPUN 
C 

RETURN 

END 

C 

c - 

c 

SUBROUTINE  USER 
C 

C  0.4  LOG  USER  CSID,  TIME  AND  DATE  TO  DIARY  FILE 
C 

INTEGER*4  INAME,  ITIME(2),  IDATE ( 3 ) 

CALL  GUINFO( 'SIGNONID' , INAME) 

CALL  TIME(4,0, ITIME) 

CALL  TIME(5,0, IDATE) 

CALL  FTNCMD( 'ASSIGN  1 8=HADI : ULOG ' , 1 9 ) 

WRITE ( 18,20)  INAME, ITIME ( 1 ) , ITIME (2 ) , ( IDATE ( J ) . d= 1 , 3 ) 
20  FORMAT ( '  MEGAERA ' , 5X , A4 , 5X , 2A4 , 5X , 3A4 ) 

RETURN 

END 

C 

- - 

c 

SUBROUTINE  LABRUN 
C 

C  1.1  LABEL  THE  RUN 
C 

C/  INSERT  COMBAS 

NAMELIST/LABELS/LABEL  1 , LABEL2 , LABEL3 . LABEL4 
C 

C  READ  IN  LABELS 

READ(NREAD, LABELS) 

C 

C  PRINT  PROGRAM  HEADING  AND  LABELS 

CALL  BLINES( 3 ) 

WRITE (NPRINT , 20) 

WRITE ( NPRINT , 2 1 ) 

WRITE ( NPRINT , 22 ) 

CALL  BLINES( 2 ) 

CALL  MESAGE ( LABEL  1 ) 

CALL  MESAGE(LABEL2) 


2.1  STEP  ON  THE  CALCULATION 

CALL  STEPON 

3 .  OUTPUT 

3.1  PERIODIC  PRODUCTION  OF  OUTPUT 
CALL  OUTPUT ( 2 ) 

4.  EPILOGUE 

4.1  TEST  FOR  COMPLETION  OF  RUN 

CALL  TESEND 

IF( .NOT .NLEND)  GO  TO  210 
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CALL  MESAGE ( LABEL3 ) 

CALL  MESAGE ( LABEL4 ) 

RETURN 

C 

C 

20  FORMAT ( 35X , ' P  ROGRAM  MEGAERA') 

21  FORMAT ( 35X , ' - ' ) 

22  FORMAT( 'O' , 34X , 'ADH' , 14X, 'AUGUST,  1980') 

C 

END 

C 

c - 

c 

SUBROUTINE  CLEAR 
C 

C  1.2  CLEAR  VARIABLES  AND  ARRAYS 
C 

c 

C/  INSERT  COMGLO 
C/  INSERT  COMELE 
C/  INSERT  COMTHM 
C/  INSERT  COMTEM 
C/  INSERT  COMCON 
C/  INSERT  COMDIM 
C 

NDX=50 
NDY  =  50 
NV=NDX 

I F ( NDY . GT . NDX )  NV=NDY 
NA=NDX*NDY 
C 

C  CLEAR  GLOBAL  VARIABLES 

I  L  =  6  *  NA  +  4  *  NV  +  7 

CALL  RESETR(ELECON.IL. 0.0) 

C 

C  CLEAR  ELECTRICAL  COMMON  BLOCKS 

I  L=  7  *  NA  +  3  *  10  +  1  *  NV 

CALL  RESETR(ECXP. IL.O.O) 

C 

C  CLEAR  THERMAL  COMMON  BLOCK 

I L=  3  *  10  +  7  *  NA 

CALL  RESETRCTHMCON. IL.O.O) 

C 

C  CLEAR  TEMPERARY  ARRAYS 

IL=  3  *  NA  +  5  *  NV 

CALL  RESETR(TEM1 . IL.O.O) 

C 

C  CLEAR  TIMESTEP  CONTROL  COMMON  BLOCK 

IL=  38 

CALL  RESETR(HTIME. IL.O.O) 

C 

C  CLEAR  INTEGER  VARIABLES 

I  L=  18 

CALL  RESET  I ( NX , I L , O ) 

C 

RETURN 

END 

C  _ 

c 
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SUBROUTINE  PRESET 
C 

C  1.3  SET  DEFAULT  VALUES 
C 

C/  INSERT  COMDIM 
C/  INSERT  COMELE 
C/  INSERT  COMCON 
C/  INSERT  COMGLO 
C 

C  THE  ELECTRICAL  AND  THERMAL  GEOMETRY  DEFAULTS  TO  INSULATORS. 

C  THIS  HAS  BEEN  PRESET  BY  SETTING  THE  ARRAYS  ELEGEO  AND 

C  THMGEO  TO  ZERO  (  AS  DONE  IN  SUBROUTINE  CLEAR  ) . 

C 

NDX =50 
NDY  =50 
EPSELE® 1 .E+6 
SCALE® 1  .0 
THICK® 1 .0 
MMAX*200 
NSPO® 1 
NERR0R®0 
NPX®20 
NPY=20 
NId=1 
NCI  =0 
NSTO® 1000 

c 

RETURN 

END 

C 

C - 

c 

SUBROUTINE  DATA 
C 

C  1.4  DEFINE  DATA  SPECIFIC  TO  RUN 
C 

C/  INSERT  COMGLO 
C/  INSERT  COMELE 
C/  INSERT  COMTHM 
C/  INSERT  COMDIM 
C/  INSERT  COMBAS 
C/  INSERT  COMCON 
C 

DATA  HHOR , HVER/ 'HORI ' , ' VERT '/. HELE . HINS/ '  ELEC  '  , 
ft  ' INSU ' / , HCTE/ ' CTEM ' / 

C 

C  NAMELISTS 

NAMELIST  /INPUT  1/ 

H  NX . NY . DELTAX , DELTAY . THICK , GEOMET , NREG 

NAMELIST  /REGION/ 

»  MINI , MAXI . MINd . MAXd . EALPHA . EBETA . THMK , 

*  THMRC , ETYPE . TTYPE , VOLTS , TEMPER 

NAMELIST  /INPUT2/ 

H  XMIN , XMAX , YMIN , YMAX . NSPO . NSTO . NPX . NPY 

NAMELIST  /INPUT3/ 

ft  TINIT.DTIME.HTIME.NTYPE.CCUR.CPOW.CTIME, 

ft  DTEMP.NId.NCI  .CVOL.EPSELE.MMAX 

C 

C  RESUMING  AN  OLD  RUN? 

IF(NLRES)  GOTO  10 
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READ  IN  FIRST  NAMELIST 
GEOMET  =HVER 

READ( NREAD , INPUT  1 , END  =900) 

NALPHA  =MAX ( NX , NY ) 

CALL  BLINES( 2 ) 

CALL  OUTR( 'THICKNESS  OF  2D  SLICE  (THICK) 
CALL  OUTH( 'GEOMETRY  OF  PROBLEM 
CALL  OUTI( 'NUMBER  OF  REGIONS  (NREG) 

I F ( NREG . GT . 9 )  GOTO  901 
NGE0=0 

IF ( GEOMET . EO.HHOR)  NGEO= 1 

READ  IN  AND  OVERLAY  DIFFERENT  REGIONS 
DO  1  K= 1 , NREG 

READ (NREAD, REGION. END=900) 

CALL  BLINES( 2 ) 

CALL  OUTI( 'REGION  NUMBER 


FIND  ELECTRICAL  AND  THERMAL 
I ET=-K 

IF(ETYPE.EO.HELE) 
I F ( ETYPE . EO.HINS) 
ITT  *-K 

IF(TTYPE. EO.HINS) 
IF(TTYPE . EQ.HCTE ) 


REGION  TYPES 


I ET  SK 
I  ET=0 


ITT =0 
ITT =K 


'  .THICK) 

' .GEOMET) 
' , NREG ) 


'  .K) 


OVERLAY  ELECTRICAL  AND  THERMAL  GEOMETRIES 
DO  2  I =MINI , MAXI 

DO  3  <J=MINJ , MAXU 

ELEGEO( I , J)=IET 
THMGEO( I , J ) =ITT 
CONTINUE 
CONTINUE 


4 


5 


C 

1 


STORE  ELECTRICAL  AND  THERMAL  PROPERTIES  OF  THE  REGION 

CALL  OUTH( 'ELECTRICAL  TYPE  '.ETYPE) 

IF(IET.GE.O)  GOTO  4 
ELEALP(K ) =EALPHA 
ELEBET(K ) =EBETA 

CALL  OUTR ( ' TEMP .  DEPENDENCE  OF  CONDUCTIVITY', 

*  EALPHA) 

CALL  OUTR( 'CONDUCTIVITY  AT  24  CELSIUS  '  . EBETA  ) 

IF(IET.GT.O)  ELEVOL ( K )= VOLTS 
IF(IET.GT.O) 

»  CALL  OUTR( 'VOLTAGE  (  BEFORE  SCALING  )  '.VOLTS) 

CALL  OUTH( 'THERMAL  REGION  TYPE  '.TTYPE) 

IF(ITT.GE.O)  GOTO  5 
THMCON( K ) =THMK 
THMCAP( K )=THMRC 

CALL  OUTR( 'THERMAL  CONDUCTIVITY 
„  THMK ) 


CALL  OUTR( 'THERMAL  HEAT  CAPACITY 

THMRC ) 

IF(ITT.GT.O)  THMTEM(K)=TEMPER 
IF(ITT.GT.O) 

CALL  OUTR ( 'CONSTANT  TEMPERATURE  '.TEMPER) 


CONTINUE 

XCOORD(  1  )  =DELTAX ( 1 )/2. 
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YC00RD( 1 ) =DELTAY ( 1 >/2 . 

DO  6  I =2. NX 

XC00RD ( I  )=XC00RD( 1-1 )+ ( DELTAX ( I -  1 )+DELTAX ( I ) )/2. 
DO  7  J=2 , NY 

YCOORD ( J ) = YC00RD( J-1 )+(DELTAY(J-1 )+DELTAY(d) )/2. 
XMIN=XC00RD (  1  ) 

XMAX=XC00RD ( NX ) 

YMIN=YCOORD(  1  ) 

YMAX= YCOORD ( NY ) 


INPUT  OUTPUT  CONTROL  VARIABLES 
RE AD ( NREAD , INPUT2 ) 

CALL  BLINES( 2 ) 

CALL  OUT  I (  ' NO .  OF  STEPS  PER  PRINTED  OUTPUT 
CALL  OUT  I ( ' NO .  OF  STEPS  PER  TAPE  STORAGE 
CALL  OUTI ( ' NO .  OF  INTERPOLATION  POINTS  IN  X 
CALL  OUT  I ( ' NO .  OF  INTERPOLATION  POINTS  IN  Y 
CALL  OUTR ( ' X  MINIMUM  FOR  INTERPOLATION 
CALL  OUTR ( ' X  MAXIMUM  FOR  INTERPOLATION 
CALL  OUTR('Y  MINIMUM  FOR  INTERPOLATION 
CALL  OUTR ( ' Y  MAXIMUM  FOR  INTERPOLATION 


' , NSPO ) 
' , NSTO ) 
'  , NPX  ) 

'  , NPY  ) 

' , XMIN ) 
' , XMAX ) 
' , YMIN ) 
'  , YMAX  ) 


INPUT  HEATING  INFORMATION 
10  READ(NREAD, INPUT3) 

IF(EPSELE.GT. 1 .E+4)  EPSELE=CVOL/ 1 . E+5 
DELT=DTEMP 

I F ( NCI . EO.O)  GOTO  902 
CALL  BLINES( 3 ) 

CALL  OUTR( ' INITIAL  TEMPERATURE  IN  CELSIUS  '.TINIT) 
CALL  OUTR( 'TOTAL  HEATING  TIME  IN  SECONDS  ' , HTIME ) 
CALL  OUTR( ' INITIAL  TIMESTEP  SIZE  IN  SECONDS DTIME ) 
CALL  OUTI( 'GRIDLINE  FOR  CURRENT  INTEGRAL  ',NCI> 

I F ( NTYPE . EO. 1 ) 

H  CALL  OUTR( 'CONSTANT  VOLTAGE  '.CVOL) 

I F ( NTYPE . E0.2) 

M  CALL  OUTR( 'CONSTANT  CURRENT  (  AMPERES  )  ' , CCUR ) 

I F ( NTYPE . EO . 3 ) 

H  CALL  OUTR( 'CONSTANT  POWER  (  WATTS  )  '.CPOW) 

CALL  OUTR( 'CHANGE  IN  TEMPERATURE  /  TIMESTEP  ' , DTEMP ) 
CALL  OUTI ( 'MAXIMUM  NO.  OF  ITERATIONS  ( MMAX ) ' , MMAX ) 
CALL  OUTR( 'ELECTRICAL  CONVERGENCE  CRITERIA  '.EPSELE) 
CALL  OUTI( 'NUMBER  OF  ITERATION  PARAMETERS  '.NALPHA) 

RETURN 


900 


901 


902 


ERROR  MESSAGES 

CALLMESAGE ( 48H  ***  NAMELIST  WAS  NOT  FOUND  DURING  DATA  INPUT  ) 

CALL  ENDRUN 
STOP 

CALLMESAGE ( 4 8H  ***  ONLY  NINE  REGIONS  MAY  BE  SPECIFIED  ) 

CALL  ENDRUN 
STOP 

CALLMESAGE ( 48H  ***  NO  SURFACE  FOR  CURRENT  INTEGRAL  IS  GIVEN  ) 

CALL  ENDRUN 
STOP 

END 


C 
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1.5  SET  AUXILIARY  VALUES 

CALL  EPARAM 
CALL  TCOEFF 

RETURN 

END 


C 

C- 

C 


C 

C 

C 

C/ 

C/ 

C/ 

C/ 

C 

C 


SUBROUTINE  INITAL 

1.6  DEFINE  PHYSICAL  INITIAL  CONDITIONS 

INSERT  COMGLO 
INSERT  COMDIM 
INSERT  COMCON 
INSERT  COMTHM 

SET  INTI AL  TEMPERATURE 
DO  1  J= 1 , NY 


DO  1 


CONTINUE 


1= 1 .NX 

TEMP( I , J)=TINIT 
K=THMGEO ( I ,U) 

IF(K.GT.O)  TEMP ( I , J ) =THMTEM( K ) 


RETURN 

END 


C 

C- 

C 


SUBROUTINE  RESUME 


C 

C  1.7  RESUME  FROM  PREVIOUS  RECORD 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMGLO 
C/  INSERT  COMELE 
C/  INSERT  COMTHM 
C/  INSERT  COMCON 
C/  INSERT  COMDIM 
C 

C  SET  DEFAULT  VALUES 

NDX=50 
NDY  =  50 
SCALE  = 1  .0 
NERR0R=0 

INPUT  RECORD  OF  RUN  GEOMETRY  AND  MATERIAL  PROPERTY  DATA 
READ ( NR E SUM . END =900 ) ELEGEO , THMGEO .DELTAX.DELTAY , THICK , 
if  XCOORD  ,  YCOORD  .  XMIN,  XMAX  ,  YMIN  ,  YMAX  .  ELEVOL  . 

M  ELEALP.ELEBET.EPSELE. THMCON , THMCAP . THMTEM , 

H  TINIT . NX , NY . NALPHA , MMAX . NREG . NGEO . NSPO . 

a  NSTO.NPX ,NPY .NIU.NCI 

READ  UNTIL  SPECIFIED  RECORD  IS  INPUT 

READ ( NRE SUM , END =900 )NSTEP , TIME , POTENT , QTHERM , TEMP , TEENG , TENG. 
M  TOENG, OENG , NR , DTI ME 

I F ( NR . LT . NREC )  GOTO  2 


C 

C 


C 

c 
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IF (NR . GT .NREC)  GOTO  900 
RETURN 
C 

C  RECORD  NOT  FOUND 

900  CALLMESAGE ( 48H  ***  SPECIFIED  RECORD  NOT  FOUND  IN  RESUME  ) 

CALL  ENDRUN 
STOP 
END 
C 

C - 

c 

SUBROUTINE  START 
C 

C  1.8  START  OR  RESTART  THE  RUN 

C 

C 

RETURN 

END 

C 

C - 

c 

SUBROUTINE  EPARAM 
C 

C  1.9  CALCULATE  ITERATION  PARAMETERS  FOR  ELECTRICAL  A.D.I.P. 

C 

C/  INSERT  COMELE 
C/  INSERT  COMDIM 
C  • 

DO  3  K=  1 , NALPHA 

ELEPAR(K)=(SIN(K*3 . 141528/(2 . *NALPHA ) ) )**2 
3  CONTINUE 

C 

RETURN 

END 

C 

- - 

c 

SUBROUTINE  TCOEFF 
C 

C  1.10  CALCULATE  THE  COEFFICENTS  FOR  THE  HEAT  EQUATION 
C 

C/  INSERT  COMGLO 
C/  INSERT  COMDIM 
C/  INSERT  COMTHM 
C 

C  COVER  THE  GRID 

DO  1  J= 1 , NY 
JP=U+1 
JM=U- 1 

DY  =DELTA Y ( U ) 

DO  2  1= 1 ,NX 
IP=I+1 
I M= I  -  1 
C 

C  CHECK  IF  THE  BLOCK  IS  IN  THE  THERMAL  DOMAIN 

K=THMGEO ( I , J ) 

IF(K.GE.O)  GOTO  2 

c 

C  FIND  THE  THERMAL  CONDUCTIVITY  AND  CLEAR  VARIABLES 

TC=THMCON( -K ) 
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DX=DELT  AX ( I ) 
TCRHS( I , J  )  =0 
TCX ( I , J ) =0 . 
TCY ( I , J ) =0 . 


FIND  COEFFICENTS  IN  THE  X  DIRECTION 
CHECK  SURROUNDING  BLOCKS 
KP=0 

IF(I.LT.NX)  KP=THMGEO( IP. J) 
KM=0 

IF(I.GT.I)  KM=THMGEO( IM.J) 
TCP=TC 

IF(KP.LT.O)  TCP=THMCON( -KP ) 
TCM=TC 

IF(KM.LT.O)  TCM=THMCON( -KM) 

RP=TC/TCP 

RM=TC/TCM 


DECIDE  WHICH  CASE  THIS  GRID  BLOCK  IS 
1.  NORMAL  BLOCK  SURROUNDED  BY  DOMAIN 

IF(KP.GE.O.OR.KM.GE.O)  GOTO  3 
H=(DX+DELTAX( IP)  )/2. 

HM=(DX+DELTAX( IM) )/2. 

TCXP ( I ,J)=(TC+TCP)/(H*( H+HM ) ) 

TCXM( I , J)=(TC+TCM)/(HM*(H+HM) ) 

CHECK  FOR  CONDUCTIVITY  DISCONTINUITIES 
IF(KP .NE .K) 

M  TCXP( I . J)=TC/( . 25* ( H+HM )*(DX+RP*DELTAX( IP ) ) ) 

'  IF(KM.NE.K) 

H  T  CXM( I , J)=TC/(  . 25* ( H+HM )*(DX+RM*DELTAX( IM) ) ) 

TCX ( I , J ) =TCXM( I , J )+TCXP ( I , d ) 

GOTO  10 


H 


H 


CONSTANT  TEMPERATURE  BOUNDARY  IN  ONE  DIRECTION 
IF(KP.EO.O.OR.KM.EO.O)  GOTO  4 
H=(DX+DELTAX( IP) )/2 . 

HM= (DX+DELTAX( IM) )/2 . 

IF(KP.GT.O)  H=DX/2 . 

IF(KM.GT.O)  HM=DX/2 . 

TCXP ( I ,J)=(TC+TCP)/(H*( H+HM ) ) 

TCXM( I , d ) = ( TC+TCM )/ ( HM* (H+HM) ) 

TCX ( I , J  )  =TCXP ( I , J )+TCXM( I . J) 

IF(KP.GT.O)  TCRHS( I . J ) =TCRHS ( I  ,d)  + 

TCXP(I.J) *THMTEM(KP ) 
IF(KP.GT.O)  TCXP ( I , J ) =0 . 

IF(KM.GT.O)  TCRHS ( I . J ) =TCRHS (  I , J  )  + 

TCXM( I , J ) *THMTEM( KM ) 


IF(KM.GT.O)  TCXM ( I , J ) =0 . 
GOTO  10 


3.  BOUNDARY  IN  ONE  DIRECTION 

IF(KP.GE.O. AND . KM . GE . 0 )  GOTO  900 
TCXM( I . J ) =0 . 

TCXP ( I . J)=0. 

IF(KP.LT.O)  TCXP( I  ,  d  )  = 

(TC+TCP  )/ (  .5*( (DX+DELTAX( IP)  )**2)  ) 
IF(KM.LT.O)  TCXM(I.d)* 

( TC+TCM )/( .  5*  ( (DX+DELTAX( IM) )**2) ) 
TCX ( I . J)=TCXM( I , d)+TCXP( I , J) 

C 
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c 

10 


CALCULATE  COEFFICENTS  IN  Y  DIRECTION 


KP  =  0 

IF(d.LT.NY) 

KM=0 

IF( J.GT.  1  ) 

TCP=TC 

I F ( KP . LT .0) 

TCM=TC 

IF(KM.LT.O) 

RP=TC/TCP 

RM=TC/TCM 


KP=THMGEO ( I , JP) 
KM=THMGEO ( I , dM) 
TCP=THMCON( -KP) 
TCM=THMCON( -KM) 


DECIDE  WHICH  CASE  THIS  GRID  BLOCK  IS 

1 .  NORMAL  BLOCK  SURROUNDED  BY  DOMAIN 

I F ( KP . GE . 0 . OR . KM . GE . O )  GOTO  5 
H=(DY+DELTAY(dP) )/2 . 

HM=(DY+DELTAY(dM) )/2 . 

TCYP( I . d)=(TC+TCP)/(H*(H+HM) ) 

TCYM( I , d )  =  ( T  C  +  TCM ) / ( HM* ( H+HM ) ) 

CHECK  FOR  CONDUCTIVITY  DISCONTINUITIES 
IF (KP . NE . K) 

H  TCYP ( I ,d)=TC/( . 25 *( H+HM ) * ( DY+DELTA Y ( UP ) *RP ) ) 

IF(KM.NE.K) 

H  TCYM( I , d)=TC/( . 25* ( H+HM )*(DY+DELTAY(dM)*RM) ) 

TCY ( I , d ) =TCYM( I , d ) +TC YP ( I , d ) 

GOTO  2 

2.  CONSTANT  TEMPERATURE  BOUNDARY 

5  IF(KP.EO.O.OR.KM.EO.O)  GOTO  6 

H=(DY+DELTAY(dP) )/2. 

HM=(DY+DELTAY(dM) )/2. 

IF(KP.GT.O)  H=DY / 2 . 

IF(KM.GT.O)  HM=DY/2 . 

TCYP( I ,d)=(TC+TCP)/ (H*( H+HM ) ) 

TCYM( I , d ) = ( TC+TCM ) / ( HM* ( H+HM ) ) 

TCY ( I , d ) =TCYM( I ,d)+TCYP( I , d ) 

IF(KP .GT .0)  TCRHS( I , d ) =TCRHS ( I  ,d)  + 

U  TCYP(I,d) *THMTEM( KP ) 

IF(KP.GT.O)  TCYP( I ,d) =0 . 

IF(KM.GT .0)  TCRHS ( I ,d)=TCRHS( I  ,  d )  + 

#  TCYM(I,d) *THMTEM( KM ) 

IF(KM.GT.O)  TCYM( I . d ) =0 . 

GOTO  2 


3.  BOUNDARY  IN  ONE  DIRECTION 

IF(KP . GE .0. AND . KM . GE .0)  GOTO  900 
TCYP ( I , d ) =0 . 

TCYM( I , d ) =0 . 

IF(KP.LT.O)  TCYP (I , d ) = 

( TC+TCP )/ ( ,5*( (DY+DELTAY (dP) )**2) ) 

IF(KM.LT.O)  TC YM ( I , d ) = 

( TC+TCM )/( . 5*( (DY+DELTAY(dM) )**2) ) 
TCY ( I , d ) =TCYM( I , d )+TCYP( I , d ) 


2  CONTINUE 

1  CONTINUE 

RETURN 

ERROR  RETURN 
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900  CALLMESAGE( 48H  ***  ERROR  RETURN  FROM  ROUTINE  TCOEF  ) 

CALL  ENDRUN 
RETURN 

END 

C 

C - - 

c 

SUBROUTINE  STEPON 
C 

C  2.1  STEP  ON  THE  CALCULATION 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMGLO 
C/  INSERT  COMCON 
C/  INSERT  COMDIM 

NSTEP=NSTEP+ 1 
C 

c 

C  IF  END  OF  HEATING  PERIOD.  RESET  ELECTRICAL  VARIABLES 

IF( ABS(TIME-HTIME ) . GT . 1 . )  GOTO  1 
NTYPE  =4 

CALL  RESETR( POTENT, NDX*NDY .0.0) 

CALL  RESETR ( QTHERM , NDX*NDY ,0.0) 

V0LTS=0 . 

P0WER=0 . 

CUR  =0. 

RES  I ST=0 . 

C 

C  IF  NO  HEATING,  SKIP  E-FIELD  AND  0  CALCULATIONS 

I  I F ( NTYPE . EO . 4 )  GOTO  14 
C 

C  CALCULATE  ELECTRICAL  CONDUCTIVITY 

CALL  ECOND 

I F ( NERROR . NE . 0 )  RETURN 
C 

C  DESCALE  POTENTIAL  AND  FIND  NEW  POTENTIAL 

S= 1 ./SCALE 

CALL  SCALER ( POTENT , NDX*NDY , S ) 

CALL  ELEPOT 

IF(NERROR .NE .0)  RETURN 
C 

C  SCALE  THE  POTENTIAL  IF  CONSTANT  CURRENT  OR  POWER 

CALL  ELECUR 
RESIST =C VOL /CUR 
GOTO  (10,11.12)  .NTYPE 

C 

C  CONSTANT  VOLTAGE 

VOLTS  =  CVOL 
POWER= VOLTS* CUR 
GOTO  13 
C 

c  CONSTANT  CURRENT 

I I  SCALE  =  CCUR /CUR 

VOLT S=CVOL* SCALE 

CALL  SCALER( POTENT, NDX*NDY. SCALE) 

POWER=CCUR* VOLTS 

CUR=CCUR 

GOTO  13 

C 

c  CONSTANT  POWER 
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12  SCALE=SQRT(CP0W/( CUR *C VOL ) ) 

V0LTS=SCALE*CV0L 

P0WER=CP0W 

CUR =CUR* SCALE 

CALL  SCALER ( POTENT , NDX*NDY , SCALE ) 

CALCULATE  THE  HEATING  RATES 

13  CALL  OCALC 

I F ( NERROR . NE . 0 )  RETURN 

FIND  THE  TIMESTEP  SIZE  FOR  THE  NEXT  HEATING  STEP 

14  DTIME=DTIME*DTEMP/DELT 

I F ( NTYPE . LE . 3 . AND . (DTIME+TIME ) . GT . (HTIME- 1 . E-2  )  ) 
H  DTIME=HTIME-TIME 

I F  (  NTYPE . EO . 4 . AND . (DTIME+TIME) . GT . ( HTIME+CTIME ) ) 
H  DTIME=HTIME+CTIME-TIME 

FIND  THE  NEW  TEMPERATURE 
CALL  TCALC 

IFINERROR .NE .0)  RETURN 

CHECK  THE  ENERGY  BALANCE 
CALL  ENGBAL 
TIME=TIME+DTIME 
RETURN 
END 


SUBROUTINE  ELEPOT 

2.2  CALCULATE  THE  ELECTRICAL  POTENTIAL  FROM  THE 
FINITE  DIFFERENCE  EQUATION,  USING  THE 
ALTERNATING  DIRECTION  IMPLICIT  PROCEDURE. 

/  INSERT  COMGLO 
/  INSERT  COMELE 
/  INSERT  COMTEM 
/  INSERT  COMDIM 

FIRST,  CALCULATE  THE  COEFFICENTS  FOR  THIS  TIME  STEP 
CALL  ELECOF 

IF(NERROR .NE .0)  RETURN 

USE  THE  POTENTIAL  OF  THE  LAST  STEP  AS  THE  FIRST  GUESS 
CALL  COP YR ( POTENT , 1 , TEM 1 , 1 , NDX*NDY ) 

RESET  THE  ITERATION  COUNTER 
MITER=0 

***  START  OF  ITERATION  LOOP  FOR  A.D.I.P.  *** 

1  MITER  =  MITER+ 1 

CHECK  IF  THE  MAXIMUM  NUMBER  OF  ITERATIONS  IS  REACHED 
IF(MITER.GT.MMAX)  GOTO  900 

CHOOSE  THE  ITERATION  PARAMETER 

L=MITER-( MITER/NALPHA ) *NALPHA 
IF(L.EQ.O)  L=NALPHA 
ALPHA=ELEPAR( L ) 


■ 
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114  1 

C 

1  142 

C 

FIRST  STAGE,  IMPLICIT  IN  X. 

1  143 

DO  2  J= 1 , NY 

1  144 

JM= J- 1 

1  145 

IF(  J.  EO.  1  )  <JM=  1 

1  146 

dP=d+1 

1  147 

IF(d.EQ.NY)  dP=NY 

1  148 

NEQU=0 

1  149 

DO  3  1=1 ,NX 

1  150 

C 

1151 

C 

IF  THE  BLOCK  IS  AN  ELECTRODE,  OR  NOT  II 

1  152 

C 

SKIP  TO  THE  NEXT  BLOCK. 

1153 

K=ELEGEO( I . J) 

1  154 

IF(K.GE.O)  GOTO  3 

1155 

C 

1156 

C 

CALCULATE  AND  ASSIGN  TRIDIAGONAL  COEFF 

1157 

NEQU=NEQU+1 

1  158 

ECSUM=EXMXP ( I . d)+EYMYP( I , d) 

1  159 

T 1 ( NEOU ) =ECXM( I , J) 

1  160 

T2(NEQU)=-(EXMXP( I , d )  +  ALPHA* 

1  161 

T3(NE0U)=ECXP( I , J) 

1  162 

T4( NEOU )=-(ECYM( I . J ) *TEM 1 ( I , 

1  163 

H  +( EYMYP( I , d)-ALPHA* 

1  164 

#  -ERHS(I.J) 

1  165 

3 

CONTINUE 

1  166 

C 

1  167 

C 

SOLVE  THE  TRIDIAGONAL  SYSTEM 

1  168 

CALL  THOMAS 

1  169 

C 

1170 

c 

PLACE  RESULTS  IN  TEM2 

1  17  1 

NEQU=0 

1172 

00  5  1=1 , NX 

1  173 

TEM2 ( I , J ) =0 . 

1  174 

K=ELEGEO( I , J) 

1  175 

IF(K.GE.O)  GOTO  5 

1  176 

NEQU=NEQU+ 1 

1177 

TEM2( I ,J)=T5(NEQU) 

1  178 

5 

CONTINUE 

1  179 

c 

1  180 

2 

CONTINUE 

118  1 

C 

1  182 

C 

START  OF  SECOND  STAGE 

1  183 

DO  6  1=1 , NX 

1  184 

IM= I -  1 

1  185 

IF( I . EQ. 1 )  I M= I 

1  186 

IP=I+1 

1  187 

IF(I.EQ.NX)  IP=I 

1  188 

NE0U=O 

1  189 

DO  7  J=1 .NY 

1  190 

C 

1  191 

C 

IF  NOT  IN  DOMAIN,  SKIP  TO  NEXT  POINT 

1  192 

K=ELEGEO( I . J) 

1  193 

IF(K.GE.O)  GOTO  7 

1  194 

c 

1  195 

c 

CALCULATE  AND  ASSIGN  TRIDIAGONAL  COEFF 

1  196 

NE0U=NE0U+1 

1  197 

ECSUM=  EXMXP ( I , J )  +  EYMYP ( I ,d) 

1  198 

T 1 ( NEOU )=ECYM(I,d) 

1  199 

T2(NE0U)=-( EYMYP( I , d)+ ALPHA* 

1200 

T3( NEOU ) =ECYP( I . d ) 

THE  DOMAIN, 


J )  *TEM 


1(1 ,dP)) 
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T4 ( NEQU )  =  -(ECXM( I , J)*TEM2( IM , J )  +  ECXP ( I , J)*TEM2( IP ,  J)  ) 

*  + ( EXMXP ( I , J) -ALPHA *ECSUM)*TEM2( I ,  d) 

*  -ERHS(I.J) 

7  CONTINUE 

C 

C  SOLVE  TRIDIAGONAL  SYSTEM 

CALL  THOMAS 
C 

C  PLACE  RESULTS  IN  TEM3 

NEQU=0 
DO  9  U= 1 , NY 

TEM3 ( I , U ) =0 . 

K=ELEGEO( I , J) 

IF(K.GT.O)  TEM3 ( I , J ) =ELEVOL ( K ) 

IF(K.GE.O)  GOTO  9 

NEQU=NEQU+1 

TEM3( I ,J)=T5( NEOU ) 

9  CONTINUE 
C 

6  CONTINUE 
C 

C  CHECK  FOR  CONVERGENCE 

CALL  MAXD I F ( 1 .EPS) 

IF(EPS.LT.EPSELE)  GOTO  10 
C 

C  NO  CONVERGENCE.  TRANSFERE  TEM3  TO  TEM 1 

CALL  COP YR ( TEM3 , 1 , TEM 1 . 1 , NDX*NDY ) 

C 

C  LOOP  BACK  TO  START  NEXT  ITERATION 

GOTO  1 
C 

C  ***  END  OF  ITERATION  LOOP  *** 

C 

C  TRANSFERE  FINAL  SOLUTION  TO  POTENT 

10  CALL  C0PYR(TEM3, 1 .POTENT, 1 ,NDX*NDY) 

RETURN 

C 

C  ERROR  MESSAGE  -  SOLUTION  FAILED  TO  CONVERGE 

900  NERROR= 1 

CALL  COPYR ( TEM3 ,  1 .POTENT.  1 , NDX *NDY ) 

RETURN 

END 

C 

C - 

c 

SUBROUTINE  ELECOF 
C 

C  2.3  CALCULATION  OF  COEFFICENTS  OF  ELECTRICAL  DIFFERENCE  EOU . 

C 

C/  INSERT  COMGLO 
C/  INSERT  COMELE 
C/  INSERT  COMDIM 
C 

C  FOR  EACH  GRID  BLOCK,  CALCULATE  THE  FOUR  COEFFICENTS 

DO  1  J= 1 , NY 

DY=DELTAY(J) 

DO  2  1=1 .NX 

C  IF  BLOCK  NOT  IN  DOMAIN,  THEN  SKIP  TO  NEXT  BLOCK 

K=ELEGEO( I . J) 

IF(K.GE.O)  GOTO  2 
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ERHS ( I , d )=0 . 

EXMXP ( I , J ) =0 . 

E YMYP( I , J ) =0 . 

DX=DELT AX ( I  ) 

S=ELEC0N( I , J) 

CALCULATE  COEFFICENTS  FOR  X  DIRECTION 
CHECK  SURROUNDING  BLOCKS 
KP  =  0 


IF  (  I  . 
KM  =  0 

LT. 

NX) 

KP=ELEGEO ( 1+1 

.J) 

I  F  (  I  . 

GT  . 

1) 

KM=ELEGEO( 1-1 

,  J ) 

SP  =  S 

I F  ( KP 
SM=S 

.  LT 

•  0) 

SP=ELECON( 1+ 

1  .d) 

IF  (KM 

.  LT 

.0) 

SM=ELECON( I - 

1  ,d) 

RP=S/SP 

RM=S/SM 

DECIDE  WHICH  CASE  THIS  GRID  BLOCK  IS 
1.  NORMAL  BLOCK  SURROUNDED  BY  DOMAIN 

IF(KP.GE.O.OR.KM.GE.O)  GOTO  3 
H= ( DX+DELTAX ( 1  +  1  )  )/2. 

HM= ( DX+DELT  AX ( I -  1 ))/2. 

ECXP( I ,d)=(S+SP)/(H*( H+HM ) ) 

ECXM( I , d)=(S+SM)/(HM*(H+HM) ) 

CHECK  FOR  CONDUCTIVITY  DISCONTINUITIES 
IF(KP.NE.K) 

H  ECXP(I , d)=S/( .25*(H+HM)*(DX+RP*DELTAX(I+1 ) ) ) 

I F ( KM . NE . K ) 

H  ECXM( I , J)=S/( .25*(H+HM)*(DX+RM*DELTAX( 1-1 ) ) ) 

EXMXP (  I  ,  d  )  =  ECXM(  I  , d)+ECXP(  I  ,  J) 

GOTO  10 


0 


ELECTRODE  IN  ONE  DIRECTION 

IF(KP.EO.O.OR.KM.EO.O)  GOTO  4 
H= ( DX+DELTAX ( 1+ 1 )  )/2. 

HM= ( DX+DELT AX ( I -  1 ) )/2. 
IF(KP.GT.O)  H=DX/2 . 

IF(KM.GT.O)  HM=DX/2 . 

ECXP(I , J)=(S+SP)/(H*(H+HM) ) 
ECXM( I ,d)=(S+SM)/( HM* ( H+HM ) ) 
EXMXP (I , d)=ECXM( I , J )+ECXP( I , J) 
IF(KP.GT.O)  ERHS ( I , d )  =  ERHS ( I , J 


IF(KP.GT.O) 
I F ( KM . GT . 0 ) 

IF(KM.GT.O) 
GOTO  10 


ECXP ( I 

ECXP ( I . J)=0. 

ERHS ( I . J)=ERHS( I 
ECXM( I 

ECXM ( I ,d)=0. 


)  + 

)  *ELEVOL ( KP ) 
)+ 

) *ELEVOL ( KM ) 


3.  BOUNDARY  IN  ONE  DIRECTION 

I F ( KP .GE .0. AND . KM . GE . 0 )  GOTO  900 
ECXM  (  I , d ) =0 . 

ECXP ( I . J)=0. 

IF(KP.LT.O)  ECXP ( I , J ) = 

(S+SP)/( .5*( (DX+DELTAX(I+1 ) )**2) ) 

IF(KM.LT.O)  ECXM( I , J ) = 

(S  +  SM)/(  .  5*  ( (DX+DELTAX( 1-1 ) ) ♦ *  2 ) ) 


' 
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EXMXP( I , d ) =ECXM( I , d)+ECXP( I  ,d) 

CALCULATE  C0EFFICENTS  IN  THE  Y  DIRECTION 
10  KP=0 

IF(d.LT.NY)  KP=ELEGE0( I ,  d+  1 ) 

KM=0 

IF(d.GT.I)  KM=ELEGE0( I , J- 1 ) 

SP  =  S 

IF(KP.LT.O)  SP=ELEC0N( I , d+  1  ) 

SM=S 

IF(KM.LT.O)  SM=ELEC0N( I , J- 1 ) 

RP=S/SP 

RM=S/SM 


DECIDE  WHI 
1 .  NORMAL 


CHECK  FOR 


H 

tt 


CH  CASE  THIS  GRID  BLOCK  IS 
BLOCK  SURROUNDED  BY  DOMAIN 
IF(KP.GE.O.OR. KM . GE . 0 )  GOTO  5 
H=(DY+DELTAY( J+1  )  )/2. 

HM=(DY+DELTAY(J-1 ) )/2. 

ECYP(I ,d)=(S+SP)/(H*( H+HM ) ) 

ECYM( I ,d)=(S+SM)/( HM* ( H+HM ) ) 

CONDUCTIVITY  DISCONTINUITIES 

I F ( KP . NE . K ) 

ECYP ( I , J ) =S/ ( . 25  * ( H+HM )*(DY+RP*DELTAY(d+1 ) ) ) 
IF(KM.NE.K) 

ECYM( I , J)=S/( . 25* ( H+HM )*(DY+RM*DELTAY(J-1 ) ) ) 
E YMYP ( I , J ) =ECYM( I ,d)+ECYP( I , J) 

GOTO  2 


2.  ELECTRODE  IN  ONE  DIRECTION 

5  IF(KP.EO.O.OR.KM.EO.O)  GOTO  6 

H=(DY+DELTAY( J+1 ) )/2 . 

HM=(DY+DELTAY( J-1 ))/2. 

IF(KP.GT.O)  H=DY/2 . 

IF(KM.GT.O)  HM=DY/2 . 

ECYP(I ,d)=(S+SP)/(H*( H+HM ) ) 

ECYM( I ,d)=(S+SM)/( HM* ( H+HM ) ) 

EYMYP ( I , J ) =ECYM( I , J )+ECYP ( I , J ) 

IF(KP.GT.O)  ERHS ( I , d ) =ERHS( I , J )+ 

H  ECYP( I ,d)*ELEVOL(KP) 

IF(KP.GT.O)  ECYP ( I , d ) =0 . 

IF(KM.GT.O)  ERHS ( I , d ) =ERHS ( I , J )+ 

H  ECYM( I ,d)*ELEVOL( KM ) 

IF(KM.GT.O)  ECYM( I , J ) =0 . 

GOTO  2 

3.  BOUNDARY  IN  ONE  DIRECTION 

6  IF(KP . GE .0. AND .KM.GE .0)  GOTO  900 

EC YM( I , d ) =0 . 

ECYP ( I , d ) =0 . 

IF(KP.LT.O)  ECYP ( I , d ) = 

„  (S+SP)/( . 5*( (DY+DELTAY(d+1 ) )**2) ) 

IF(KM.LT.O)  ECYM( I , d ) = 

H  ( S+SM  )/(  . 5* ( ( DY+DE LTAY ( J- 1 ) )**2) ) 

EYMYP( I ,d)  =  ECYM( I ,d)  +  ECYP(  I ,d) 


2  CONTINUE 

1  CONTINUE 

RETURN 
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c 

c 

900 


C 

C 

C 


ERROR  RETURN 
NERR0R=4 
RETURN 
END 


SUBROUTINE  THOMAS 
C 

C  2.4  SOLVE  A  TRIDIAGONAL  SYSTEM  BY 
C 

C/  INSERT  COMTEM 
C/  INSERT  COMDIM 
C 


THOMAS'S  ALGORITHM 


DIMENSION  W( 50) , G( 50) 

W ( 1)=T3(1)/T2(1) 

G( 1 )=T4( 1  )/T2 (  1  ) 


DO  1  1=2, NEOU 
IM= 1-1 

DEN0M=T2( I)-T1(I)*W(IM) 

W( I ) =T3 ( I )/DENOM 

G( I )=(T4( I )-T1( I )*G( IM) )/DENOM 

1  CONTINUE 
C 

T5( NEOU ) =G( NEOU ) 

DO  2  I =2, NEOU 

I  1 =NEOU+ 1-1 
IP=I 1+1 

T5( I  1 )=G( I  1 )-W( I  1 )*T5( IP) 

2  CONTINUE 
C 

RETURN 

END 

C 

c - - 

c 

SUBROUTINE  MAXDI F ( ICODE . EPS ) 

C 

C  2.5  FIND  THE  MAXIMUM  DIFFERENCE  BETWEEN  TEM3 ,  TEM 1 
C 

C/  INSERT  COMTEM 
C/  INSERT  COMDIM 
C/  INSERT  COMGLO 
C 
C 

EPS=0. 

DO  1  J=1,NY 

DO  2  1=1 , NX 

I  F ( ICODE . EO.  1 )  K  =  ELEGEO( I , J ) 

IF( ICODE .E0.2)  K=THMGEO( I , J) 
IF(K.GE.O)  GOTO  2 
DI F  = ABS ( TEM3 ( I . J)-TEM1(I , J) ) 

I  F (  ICODE . EO . 2 ) 

H  DIF=DIF/ABS(TEM1  (  I  ,«J)  ) 

EPS=AMAX 1 (DIF , EPS ) 

2  CONTINUE 

1  CONTINUE 
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RETURN 

END 

C 

c - 

c 

SUBROUTINE  ECOND 
C 

C  2.6  CALCULATE  THE  ELECTRICAL  CONDUCTIVITY 
C 

C  THE  ELECTRICAL  CONDUCTIVITY  IS  ASSUMED  TO  BE  TEMPERATURE 

C  DEPENDENT.  ACCORDING  TO  THE  FORMULA; 

C 

C  SIGMA (T)=BETA( 1 + ALPHA ( T-24 ) ) 

C 

C  WHERE  T  IS  IN  DEGREES  CELCIUS. 

C 

C/  INSERT  COMGLO 
C/  INSERT  COMELE 
C/  INSERT  COMDIM 
C 

c 

DO  1  U= 1 , NY 

DO  2  1=1 .NX 

K=ELEGEO( I , J ) 

IF(K.GE.O)  GOTO  2 
C 

C  FIND  THE  COEFFICENTS  FOR  THE  CONDUCTIVITY 

K  =  -K 

ALPHA=ELEALP( K ) 

BETA=ELEBET(K) 

T  =TEMP ( I ,d) 

C 

ELECON(I . J)=BETA*( 1 . +ALPHA* ( T-24 . )) 

C 

2  CONTINUE 

1  CONTINUE 
RETURN 
END 
C 

C - 

c 

SUBROUTINE  QCALC 
C 

C  2.7  CALCULATE  THE  HEATING  RATE  AT  EACH  GRIDPOINT. 

C 

C/  INSERT  COMDIM 
C/  INSERT  COMGLO 
C 

C  FOR  EACH  GRID  POINT 

DO  80  d=1 .NY 
dP=d+1 
JM=U- 1 

DO  81  1=1 .NX 
I M  =  I  —  1 
IP=I+1 

QTHERM ( I , J ) =0 . 

C 

C  CHECK  IF  THIS  GRIDBLOCK  IS  IN  THE  DOMAIN 

K=ELEGEO( I .d) 

IF(K.GE.O)  GOTO  81 
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C 

C  CHECK  BLOCKS  IN  THE  X  DIRECTION 

KP  =  0 

IF(I.LT.NX)  KP  =  ELEGEO(IP,vJ) 

KM  =  0 

IF(I.GT.I)  KM=ELEGEO(IM, J) 

C 

C  CALCULATION  WHEN  SURROUNDING  BLOCKS  ARE  DOMAIN  BLOCKS 

I F ( KP . EO . 0 . OR . KM . EO . 0 . OR . ( KP . NE . K . AND . KM . EQ . K . AND 
H  .KP.LT.O).OR.( KM . NE . K . AND . KP . EO . K . AND . KM . LT . O) )GOTO  2 

H= ( DELTAX ( I )+DELTAX( IP) )/2. 

HM= ( DELT  AX ( I )+DELTAX(IM))/2. 

IF(KP.GT.O)  H=DELTAX( I )/2 . 

IF(KM.GT.O)  HM=DELTAX( I )/2 . 

HS=H**2 
HMS  =  HM*  *2 

EX=  (  POTENT  ( I P  ,  J  )  *HMS+  ( HS-HMS  )  ’"POTENT  ( I  ,  J  ) 

*  -POTENT  (IM,J)*HS)/(  H*HM*  ( H+HM )  ) 

GOTO  3 

C 

C  SPECIAL  CASES 

2  IF(KP.GE.O. AND . KM . GE . 0 )  GOTO  900 

IF((KP.EQ.O. AND . KM . LT . 0 ) . OR . KM . EQ . K ) 

H  EX= ( 2. /( DELTAX ( I  )+DELT AX ( I M ) ) ) * ( POTENT ( I , J ) -POTENT ( IM, J) ) 

I F ( ( KP . LT . 0 . AND . KM . EO . 0 ) . OR . KP . EO . K ) 

*  EX  =  (2./(DELTAX( I  )+DELTAX( IP) ) )*( POTENT ( IP , J ) -POTENT ( I , d ) ) 

C 

3  CONTINUE 
C 

C  CHECK  BLOCKS  IN  THE  Y  DIRECTION 

KP=0 

IF(J.LT.NY)  KP=ELEGEO( I ,JP) 

KM  =  0 

IF(J.GT.I)  KM=ELEGEO( I , JM ) 

C 

C  CALCULATION  WHEN  SURROUNDED  BY  DOMAIN  BLOCKS 

I F ( KP . EO . 0 . OR . KM . EO . 0 . OR . ( KP . NE . K . AND . KM . EO . K . AND 
M  . KP . LT . 0 ) . OR . ( KM . NE . K . AND . KP . EQ . K . AND . KM . LT . 0 ) )GOTO  4 

H=(DELTAY(J)+DELTAY(JP) )/2. 

HM= ( DELTA Y( J)+DELTAY( JM) )/2 . 

IF(KP.GT.O)  H=DELTAY( J)/2. 

IF(KM.GT.O)  HM=DELTA Y ( J ) /2 . 

HS=H**2 

HMS=HM**2 

EY= ( POTENT ( I , JP ) *HMS+ ( HS-HMS ) * POT ENT ( I . J) 

H  -POTENT ( I , UM ) *HS ) / ( H*HM* ( H+HM ) ) 

GOTO  5 
C 

C  SPECIAL  CASES 

4  IF(KP.GE.O.AND.KM.GE.O)  GOTO  900 

I F ( ( KP . EO . 0 . AND . KM . LT . 0 ) .OR.KM.EQ.K) 

H  EY=(2./(DELTAY(J)+DELTAY(JM) ) )*( POTENT ( I , J ) -POTENT ( I . JM) ) 

I F ( ( KP . LT . 0 . AND . KM . EO . 0 ) . OR . KP . EO . K ) 

»  EY=(2./(DELTAY(J)+DELTAY(JP) ) )  * ( POTENT ( I ,JP) -POTENT ( I , J ) ) 

C 

5  CONTINUE 
C 

C  CALCULATE  HEATING 

OTHERM( I , J)=ELECON( I ,J)*(EX**2+EY**2) 

C 
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81  CONTINUE 

80  CONTINUE 

RETURN 
C 

C  ERROR  RETURN 

900  NERR0R=2 

RETURN 

END 

C 

C - 

c 

SUBROUTINE  TCALC 


C 

C 

C 

C/ 

C/ 

C/ 

C/ 

C/ 

C 

C 

C 


8 

C 

C 


C 

C 


c 

c 


2 


2.8  SOLVE  THE  TEMPERATURE  EQUATION 

INSERT  COMGLO 
INSERT  COMTHM 
INSERT  COMTEM 
INSERT  COMCON 
INSERT  COMDIM 

COPY  THE  TEMPERATURE  INTO  TEM 1 

CALL  COPYR ( TEMP , 1 , TEM 1 , 1 , NDX*NDY ) 


M  =  0 
M  =  M+  1 


START  OF  FIRST  STAGE 
DO  1  J= 1 , NY 
dM  =  d~ 1 

I F ( J. EQ. 1 )  JM= 1 
JP=d+1 

IF(J.EQ.NY)  JP=NY 

NEQU=0 

DO  2  1=1 .NX 

IF  BLOCK  IS  NOT  IN  DOMAIN.  SKIP  TO  NEXT  BLOCK 
K=THMGEO ( I . J) 

IF(K.GE.O)  GOTO  2 
F1=THMCAP(-K)/(DTIME/2. ) 

ASSIGN  TRIDIAGONAL  COEFFICENTS 
NEQU=NEQU+1 
T 1 (NEQU)=TCXM( I . J) 

T2 ( NEQU )=-(TCX( I ,d)  +  F1  ) 

T3(NEQU)=TCXP( I ,d) 

T4(NEQU)-(TCYM(I.J)*TEM1(I.dM)+TCYP(I.J)*TEM1(I.JP)) 
,  + ( TCY ( I ,J)-F 1 ) *TEM 1 ( I , J ) -TCRHS ( I .d) 

,  -QTHERM(I.d) 

CONTINUE 


C 

C 


SOLVE  TRIDIAGONAL  SYSTEM 
CALL  THOMAS 


C 

c  PLACE  RESULTS  IN  TEM2 

NEQU=0 
DO  3  1=1 .NX 

TEM2 ( I , d ) =0 . 

K  =  THMGEO( I . d  ) 
IF(K.GE.O)  GOTO  3 


' 

—  X 


. 


1621 

1622 

1623 

1624 

1625 

1626 

1627 

1628 

1629 

1630 

1631 

1632 

1633 

1634 

1635 

1636 

1637 

1638 

1639 

1640 

1641 

1642 

1643 

1644 

1645 

1646 

1647 

1648 

1649 

1650 

1651 

1652 

1653 

1654 

1655 

1656 

1657 

1658 

1659 

1660 

1661 

1662 

1663 

1664 

1665 

1666 

1667 

1668 

1669 

1670 

1671 

1672 

1673 

1674 

1675 

1676 

1677 

1678 

1679 

1680 


no  ooooooo  noon  noon  no 


97 


NEQU=NEQU+1 

TEM2 ( I , J ) =T5( NEQU ) 

3  CONTINUE 
C 

1  CONTINUE 

START  OF  SECOND  STAGE 
DO  4  1=1 ,NX 
IP=I+1 

IF(I.EQ.NX)  I P=NX 
IM= 1-1 

IF(I.EQ.I)  IM=I 

NEQU=0 

DO  5  J= 1 , NY 

IF  BLOCK  NOT  IN  DOMAIN,  SKIP  TO  NEXT  BLOCK 
K=THMGEO ( I , J ) 

IF(K.GE.O)  GOTO  5 

ASSIGN  TRIDIAGONAL  COEFFICENTS 

F1=THMCAP(-K)/(DTIME/2. ) 

NEQU=NEQU+1 
T 1 ( NEOU ) =TCYM( I , J ) 

T2 ( NEQU )=-(TCY( I , J)+F1 ) 

T3 ( NEOU )=TCYP(I,J) 

T4 ( NEOU )=-(TCXM( I , J  )  *TEM2 ( IM,J)+TCXP( I , J )*TEM2( IP , J) ) 

#  + ( TCX ( I  ,  d  )  -  F  1  ) *TEM2 ( I , J ) -TCRHS ( I , J ) 

#  -OTHERMU.U) 

5  CONTINUE 

SOLVE  THE  TRIDIAGONAL  SYSTEM 
CALL  THOMAS 

PLACE  THE  RESULTS  IN  TEM3 
NEQU=0 
DO  6  J=1 ,NY 

TEM3 ( I , J ) =0 . 

K=THMGEO ( I , J ) 

IF(K.GT.O)  TEM3  (  I  ,  J )  =  THMTEM(  K  ) 

IF(K.GE.O)  GOTO  6 
NEQU  =  NEQU+  1 
TEM3( I . J)=T5(NEQU) 

6  CONTINUE 

4  CONTINUE 

FIND  THE  MAXIMUM  RELATIVE  CHANGE  IN  TEMPERATURE 
CALL  MAXDIF( 2 .DELT) 

IF  CHANGE  LESS  THEN  TOLERANCE  EXIT  OUT 
IF(DELT.LE.DTEMP)  GOTO  7 

CHANGE  TOO  BIG,  CHANGE  DTIME 
I F ( M . GT . 5 )  GOTO  900 

DTIME=DTIME*DTEMP/DELT 
DTIME=DTIME-DTIME* .01 
GOTO  8 

COPY  NEW  TEMPERATURE  TO  TEMP 

7  CALL  C0PYR(TEM3, 1 .TEMP, 1 ,NDX*NDY) 
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c 

c 

c 

900 


C 

C 

c 


RETURN 

ERROR  RETURN 
NERR0R=3 
RETURN 
END 


C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c/ 

c/ 

c/ 

c/ 

c 

c 


SUBROUTINE  ENGBAL 

2.9  CALCULATE  THE  ENERGY  BALANCE  FOR  THE  TIMESTEP 

THREE  METHODS  ARE  USED  TO  CALCULATE  THE  ENERGY  INPUT 

1.  ELECTRICAL  POWER  *  TIME 

2.  HEATING  RATE  *  VOLUME  *  TIME  -  HEAT  FLOW  OUT 

3.  TEMPERATURE  RISE  *  VOLUME  *  HEAT  CAPACITY 

THE  LAST  TWO  METHODS  ARE  CALCULATED  SEPERATELY  FOR  EACH  REGION. 

INSERT  COMGLO 
INSERT  COMTHM 
INSERT  COMCON 
INSERT  COMDIM 

CALCULATE  ELECTRICAL  ENERGY 
DE ENG  =  POWER  *  DTIME 
TEENG=TEENG  +  DEENG 


C 

C  CALCULATE  HEATING  ENERGY  AND  TEMPERATURE  RISE  ENERGY 

DQENG=0 . 

TTENG=0 . 

DO  111  K=1 , 10 

111  TENG( K ) =0 . 

DO  1  J= 1 , NY 

DY  =DELTAY ( J ) 

DO  2  1=1, NX 

V=DY*DELTAX( I )*THICK 


C 

C  CHECK  IF  BLOCK  IN  TEMPERATURE  DOMAIN 

K=THMGEO ( I . J) 

IF(K.GE.O)  GOTO.  2 
C 

C  CALCULATE  HEAT  PRODUCED  MINUS  HEAT  FLOW  OUT  TO  CONSTANT 

C  TEMPERATURE  BOUNDARIES. 

DQ=QTHERM( I , J)*V*DTIME 
I F ( I . EO . NX )  GOTO  3 
KP=THMGEO( 1+ 1 . J) 

IF(KP.LE.O)  GOTO  3 

D0 1  =  (TEMP ( I , J ) -THMTEM( KP ) ) *THMCON( K)*DTIME* 
H  2 . *DY*THICK/DELTAX( I ) 

DQ=DQ-DQ1 

3  IF(I.EO.I)  GOTO  4 

KM=THMGEO ( 1-1 . J) 

IF(KM.LE.O)  GOTO  4 

DO  1  =  ( TEMP ( I , J ) -THMTEM( KM ) ) *THMCON( K)* DTIME* 
2 . *DY*THICK/DELTAX( I ) 

DQ=DQ-DQ1 

IF(J.EQ.NY)  GOTO  5 
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KP=THMGE0( I , J+1 ) 

IF(KP.LE.O)  GOTO  5 

D0 1  *  (TEMP ( I .d)-THMTEM(KP) ) *THMCON( K ) *DTIME * 

#  2 . *DELTAX ( I ) *THICK/DY 

DQ=DQ-DQ1 

5  IF(J. EQ. 1 )  GOTO  6 

KM=THMGEO( I , J- 1 ) 

IF(KM.LE.O)  GOTO  6 

D0 1  =  (TEMP ( I , J)-THMTEM(KM) )*THMCON(K)*DTIME* 

*  2.*DELTAX(I) *THICK/DY 
DQ=DQ-DQ1 

6  CONTINUE 
DQENG=DQENG+DQ 
OENG( -K ) =QENG( -K )+DQ 

C 

C  CALCULATE  ENERGY  INDICATED  BY  TEMPERATURE  RISE. 

DT=( TEMP( I . J)-TINIT)*V*THMCAP(-K) 

TTENG=TTENG+DT 
1 ENG( -K ) =TENG( -K )+DT 
C 

2  CONTINUE 

1  CONTINUE 

TOENG=TOENG+DOENG 

RETURN 

END 

C 

C - 

c 

SUBROUTINE  ELECUR 
C 

C  2.10  FIND  THE  CURRENT  THROUGH  A  SURFACE  PERPENDICULAR  TO 
C  THE  X  OR  Y  DIRECTION 

C 

C  IF  NI J= 1 ,  THE  SURFACE  IS  PERPENDICULAR  TO  X,  ALONG  THE 

C  GRIDLINE  I =NCI . 

C  OTHERWISE,  THE  INTEGRAL  IS  PERPENDICULAR  TO  Y,  ALONG  THE 

C  GRIDLINE  U=NCI . 

C 

C/  INSERT  COMGLO 
C/  INSERT  COMDIM 
C/  INSERT  COMCON 
C 

C  CHECK  ORIENTATION  OF  INTEGRAL. 

I F ( NI J . NE . 1 )  GOTO  10 
C 

C  PERPENDICULAR  TO  X 

H= ( DELTAX ( NC I )+DELTAX(NCI+ 1 ) )/2 . 

HM= (DELTAX(NCI )+DELTAX ( NCI  -  1 ) )/2 . 

HS  =  H*  *2 
HMS=HM**2 
CUR=0 . 

DO  1  J= 1 , NY 

S=ELECON( NCI . J) 

E  =  ( HMS  *POTENT ( NC I +  1 , J )  +  ( HS-HMS ) *POTENT ( NCI . J ) 
H  -  HS*P0TENT(NCI-1 .J) ) / ( H*HM* ( H+HM ) ) 

CUR=CUR+E*S*DELTAY ( U ) 

1  CONTINUE 

CUR=ABS(CUR*THICK) 

RETURN 

C 
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C  PERPENDICULAR  TO  Y 

10  H= (DELTAY ( NCI )+DELTAY ( NCI + 1 ) )/2. 

HM= (DELTAY( NCI )+DELTAY ( NCI -  1 ))/2. 

HS  =  H*  *2 
HMS=HM**2 
CUR=0 . 

DO  2  1=1 , NX 

S=ELECON( I , NCI ) 

E=(  HMS*POTENT ( I , NCI+ 1 )  +  ( HS-HMS ) *POTENT ( I , NCI ) 

*  -  HS*P0TENT(I.NCI-1))/(H*HM*(H+HM)) 

CUR=CUR+E  *5*DELT  AX ( I ) 

2  CONTINUE 

CUR=ABS ( CUR*THI CK ) 

RETURN 

END 

C 

- - — - 

SUBROUTINE  OUTPUT( ICODE ) 

C 

C  3.1  CONTROL  THE  OUTPUT 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMDIM 
C/  INSERT  COMGLO 
C/  INSERT  COMCON 
C 

C  BRANCH  TO  INITIAL,  MIDDLE,  AND  FINAL  OUPUT  ROUTINES 

GOTO  (  1 ,2,3) , ICODE 
C 

C  INITIAL  OUTPUT 

1  IF(NLRES)  RETURN 
CALL  PAGE 

C  OUTPUT  ELECTRICAL  GRID  STRUCTURE 

CALL  OUTGRD ( 1 ) 

CALL  PAGE 

C  OUTPUT  THERMAL  GRID  STRUCTURE 

CALL  OUTGRD( 2 ) 

C  IF  TAPE  STORAGE  IS  REQUESTED.  STORE  GRID  ETC.  ON 

IF(NSTO.LT. 200 )  CALL  OUTTAP( 1 ) 

RETURN 
C 
C 

C  PERIODIC  OUTPUT 

C  CHECK  IF  PRINTED  OUTPUT  THIS  TIMESTEP 

2  I F ( MOD( NSTEP.NSPO) .NE .0. AND .NSTEP .NE . 1 )  GOTO 

4  CALL  PAGE 

CALL  OUTI ( 'TIMESTEP  NUMBER 
IF(NSTO.LE. 200 ) 

A  CALL  OUTI  (  'CURRENT  RECORD  NUMBER 

CALL  BLINES( 1 ) 

CALL  OUTR ( 'TIME  AT  END  OF  STEP 
CALL  OUTR ( 'SIZE  OF  TIMESTEP 
CALL  BLINES( 1 ) 

I F ( NTYPE . EQ . 4 )  GOTO  6 

CALL  OUTI ('NO.  OF  ITERATIONS  OF  A.D.I.P. 

CALL  OUTR( 'VOLTAGE  BETWEEN  ELECTRODES 
CALL  OUTR ( ' POWER  INPUT  TO  DOMAIN 
CALL  OUTR( 'CURRENT  BETWEEN  ELECTRODES 
CALL  OUTR( 'RESISTANCE  BETWEEN  ELECTRODES 
CALL  BLINES( 1 ) 


TAPE 


' .NSTEP) 

' , NREC ) 

' .TIME ) 

' .DTIME ) 


' .MITER) 

' .VOLTS) 

' , POWER ) 

' .CUR) 

' .RESIST) 


■ 

' 


1861 

1862 

1863 

1864 

1865 

1866 

1867 

1868 

1869 

1870 

1871 

1872 

1873 

1874 

1875 

1876 

1877 

1878 

1879 

1880 

1881 

1882 

1883 

1884 

1885 

1886 

1887 

1888 

1889 

1890 

1891 

1892 

1893 

1894 

1895 

1896 

1897 

1898 

1899 

1900 

1901 

1902 

1903 

1904 

1905 

1906 

1907 

1908 

1909 

1910 

1911 

1912 

1913 

1914 

1915 

1916 

1917 

1918 

1919 

1920 
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6  WRITE ( NPRINT , 20) 

WRITE ( NPRINT , 22 )  DEENG , DQENG 
WRITE ( NPRINT , 2 1 )  TEENG , TOENG , TTENG 
DO  100  K= 1 , NREG 

WRITE (NPRINT, 23)  K , QENG( K ) , TENG( K ) 
lOO  CONTINUE 

CALL  BLINES ( 1 ) 

I F ( NTYPE . EQ . 4 )  GOTO  7 
CALL  OUTINT ( 4 , 1 ) 

I F ( ( 25+NPY ) . GT . 45 )  CALL  PAGE 
CALL  OUTINT ( 5 , 1 ) 

I F ( ( 25+2 *NPY ) . GT . 50)  CALL  PAGE 

7  CALL  OUTINT ( 3 , 1 ) 

C 

C  PERIODIC  TAPE  OUTPUT 

5  I F ( NSTO . LE . 200 . AND . ( MOD ( NSTEP . NSTO ) . EO . 0 . OR . ICODE . EO . 3 ) ) 

#  CALL  OUTTAP( 2 ) 

RETURN 

C 

C 

20  F0RMAT(4X, 'ENERGY  BALANCE  (  JOULES  )') 

21  F0RMAT(4X, 'CUMULATIVE :  ELECTRICAL  ENERGY  =',1PE12.4, 

H  5X ,  ' HEAT  GENERATED  =  ' .  1  PE  1 2 . 4 , 5X . 

#  'TEMP.  RISE  ENERGY  =',1PE12.4) 

22  F0RMAT(4X,  'TIMESTEP  :  ELECTRICAL  ENERGY  =  ' , 1  PE  1 2 . 4 , 5X , 

#  'HEAT  PRODUCED  =',1PE12.4) 

23  F0RMAT(8X .  'REGION  = ' , 1 3 . 5X ,  ' HE AT  GENERATED  = '  ,  1  PE  1 2 . 4 , 5X , 

#  'TEMP.  RISE  ENERGY  =',1PE12.4) 

C 

C  FINAL  OUTPUT 

3  IF(MOD(NSTEP .NSPO) . NE .0. AND .NSTEP . NE . 1 )  GOTO  4 

I F ( MOD (NSTEP , NSTO ) . NE . 0 )  GOTO  5 
RETURN 
END 


C 

C- 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c/ 

c/ 

c/ 

c 

c 


SUBROUTINE  OUTGRD( ICODE ) 

3.2  OUTPUT  ONE  OF  THE  GLOBAL  VARIABLE  ARRAYS 


ICODE  INDICATES  WHICH  ARRAY  IS  TO  BE  OUTPUT 

1  -  ELEGEO,  THE  ELECTRICAL  GEOMETRY 

2  -  THMGEO ,  THE  THERMAL  GEOMETRY 

3  -  TEMP  .  THE  TEMPERATURE 

4  -  POTENT,  THE  ELECTRICAL  POTENTIAL 

5  -  QTHERM,  THE  HEATING 

6  -  ELECON,  THE  ELECTRICAL  CONDUCTIVITY 


INSERT  COMBAS 
INSERT  COMDIM 
INSERT  COMGLO 


PRINT  OUT  HEADING 
I F ( I  CODE . EO .  1 ) 
I F ( I  CODE . EO . 2 ) 
IF( ICODE . EO. 3) 
I F ( I CODE . EQ • 4  ) 
IF( ICODE .EQ.5) 
I F ( I  CODE . EO . 6 ) 


WR ITE ( NPR I  NT , 20) 
WRI TE ( NPR I  NT , 2 1 ) 
WR I TE ( NPR I  NT , 22 ) 
WRITE ( NPRINT , 23 ) 
WRITE ( NPRINT , 24 ) 
WR I TE ( NPR I  NT , 25 ) 


■ 


1921 

1922 

1923 

1924 

1925 

1926 

1927 

1928 

1929 

1930 

1931 

1932 

1933 

1934 

1935 

1936 

1937 

1938 

1939 

1940 

1941 

1942 

1943 

1944 

1945 

1946 

1947 

1948 

1949 

1950 

1951 

1952 

1953 

1954 

1955 

1956 

1957 

1958 

1959 

1960 

1961 

1962 

1963 

1964 

1965 

1966 

1967 

1968 

1969 

1970 

197  1 

1972 

1973 

1974 

1975 

1976 

1977 

1978 

1979 

1980 


102 


c 

c 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 
C 


C 


C 


3 

2 

1 


C 

C 

4 


5 


C 

c — 
c 


IF( ICODE.LT. 3)  GOTO  4 
FORMAT  STATEMENTS 

FORMAT( '0' , 10X, 'ELECTRICAL  REGIONS  OF  THE  GRID  (NOT  TO  SCALE)') 
FORMAT( '0' . 10X, 'THERMAL  REGIONS  OF  THE  GRID  (NOT  TO  SCALE)') 
FORMAT( '0' , 10X. 'TEMPERATURE  IN  DEGREES  CELSIUS') 

FORMAT( '0' . 10X, 'ELECTRICAL  POTENTIAL  IN  VOLTS') 

FORMAT( '0' . 10X, 'HEATING  RATE  IN  J0ULES/SEC0ND-METER**3 ' ) 

FORMAT( '0' . 10X, 'ELECTRICAL  CONDUCTIVITY  IN  MHOS/METER') 

FORMAT ( ' O  COLUMN  '.10(8X,I2)) 

FORMAT ( '  X  ' , 10(3X,F7 .3) ) 

FORMAT ( '  ROW  Y  ') 

FORMAT ( '  ' ,13, 1X.F7.3, 10(3X,F7.2)) 

FORMAT ( '  '  ,13,  1X.F7.3,  10(3X,F7.0)) 

FORMAT ( '  ' ,13, 1X.F7.3, 10(3X,F7.4)) 

FORMAT ( '0' ) 

FORMAT ( '  ' ,12. 5X, 50(12)) 

FORMAT( 'OCOLUMN  ',50(12)) 

FORMAT ( 'OROW' ) 

DO  1  I  1  =  1 , NX,  10 

I T  =MI N0( I  1+9, NX) 

DO  2  J2= 1 , NY , 50 

JI =MINO( U2+49.NY ) 

PRINT  PAGE  HEADING 

WRITE(NPRINT,26)  (I, 1  =  11, IT) 

WRITE ( NPRINT , 27 )  (XCOORD( I ) . 1=11 , IT) 

WRITE ( NPR INT , 28 ) 

PRINT  EACH  ROW 

DO  3  JM= J2 , JI 
U  =  NY+ 1  - JM 

I F ( I  CODE . EO . 3  ) 

ft  WRITE  (  NPRINT  ,  29  )  U  ,  YCOORD(  J )  ,  ( TEMP ( I  ,  d  )  ,  I  =  I  1  ,  IT  ) 

I F ( ICODE . EO . 4  ) 

tt  WRITE  ( NPRINT  ,  29  )  J  ,  YCOORD(  J  )  ,  ( POTENT  ( I  ,  U  )  ,  I  =  I  1  ,  IT  ) 

I F ( ICODE . EO . 5 ) 

ft  WRITE  (  NPRINT  ,  30)  U  ,  YCOORD  (  J )  ,  ( QTHERM(  I, J), 1  =  11,  IT) 

I F ( I  CODE . EO • 6 ) 

M  WR I TE ( NPRINT , 3 1 )  J , YCOORD ( J ) , ( ELECON( I , J ) , I  =  I  1 . IT ) 

CONTINUE 

WRITE ( NPRINT , 32 ) 

CONTINUE 

CONTINUE 

RETURN 

OUTPUT  THE  GEOMETRY  DEFINING  ARRAYS 
WRITE ( NPRINT , 35 ) 

DO  5  J=  1 , NY 

JM=NY+ 1-J 
I F ( ICODE . EO. 1  ) 

#  WRITE (NPRINT , 33 )  JM , ( ELEGEO( I , JM ) , I = 1 , NX ) 

I F ( ICODE . E0.2) 

#  WRI TE ( NPRINT , 33 )  JM . ( THMGEO( I , JM) , I = 1 , NX ) 

CONTINUE 

WR I TE ( NPR I  NT . 34 )  ( I . I  =  1 , NX ) 

RETURN 

END 


. 

x 


1981 

1982 

1983 

1984 

1985 

1986 

1987 

1988 

1989 

1990 

1991 

1992 

1993 

1994 

1995 

1996 

1997 

1998 

1999 

2000 

2001 

2002 

2003 

2004 

2005 

2006 

2007 

2008 

2009 

2010 

201  1 

2012 

2013 

2014 

2015 

2016 

2017 

2018 

2019 

2020 

2021 

2022 

2023 

2024 

2025 

2026 

2027 

2028 

2029 

2030 

2031 

2032 

2033 

2034 

2035 

2036 

2037 

2038 

2039 

2040 


non  non  o  non  o  non 
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FUNCTION  ALINTP(X, Y.Z.Z1 .Z2.M.N) 

C 

C  3.3  INTERPOLATE  VALUES  FROM  GRID 
C 

C  THIS  FUNCTION  INTERPOLATES  A  VALUE  AT  THE  POINT  (X,Y)  FROM 

C  THE  MATRIX  Z.  THE  INTERPOLATION  IS  LINEAR. 

C 

DIMENSION  Z ( 50 , 50 ) , Z 1 ( 50) , Z2 ( 50 ) 

FIND  FIRST  X  CO-ORDINATE  NOT  LESS  THAN  X 

DO  1  11=2, M 
i  =  i  i  - 1 

I F ( X-Z  1  ( 1 1  )  )  3 ,  1 ,  1 


1  CONTINUE 
3  XL=X-Z1(I) 

XR  =  Z  1  ( 1+1 )-X 
XBASE  =XR  +  XL 

FIND  FIRST  Y  CO-ORDINATE  NOT  LESS  THAN  Y 

DO  2  Jd=2,N 
U= J J- 1 

IF(Y-Z2(dd)  )4,2,2 


2  CONTINUE 
4  YL=Y-Z2(J) 

YR  =  Z2( J+1  )-Y 
YBASE=YR+YL 

BILINEAR  INTERPOLATION 

ALINTP= (YL*(XR*Z( I , J+ 1 )+XL*Z( 1+1 , J+1 )) 

M  +  YR* ( XR*Z( I ,d)+XL*Z(I+1 , J) ) )/(XBASE*YBASE) 

RETURN 
END 


C 

C 

C 

C 

C 

C 

C 

c/ 

c/ 

c/ 

c/ 

c/ 

c/ 

c 

c 


SUBROUTINE  OUTINT( ICODE , IP ) 

3.4  INTERPOLATE  ONE  OF  THE  GLOBAL  VARIABLES 

ICODE  -  DETERMINES  WHICH  VARIABLE  IS  INTERPOLATED 
IP  =  1,  IF  PRINTING  DESIRED. 

RESULTS  OF  I TERPOL AT  I  ON  STORED  IN  TEM3 


INSERT  COMGLO 
INSERT  COMTEM 
INSERT  COMDIM 
INSERT  COMBAS 
INSERT  COMELE 
INSERT  COMTHM 


TX  =XMAX -XMI N 
TY=YMAX-YMIN 
SX  =TX/ ( NPX -  1  .  ) 
SY=TY/(NPY-1 . ) 


SEE  OUTGRD  ) 
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204  1 

C 

2042 

C 

204  3 

2044 

2045 

2046 

2047 

C 

2048 

C 

2049 

C 

2050 

2051 

2052 

2053 

2054 

2055 

2056 

2057 

2 

2058 

1 

2059 

C 

2060 

C 

2061 

2062 

c 

2063 

c 

2064 

c 

2065 

2066 

2067 

2068 
2069 

c 

2070 

c 

2071 

22 

2072 

23 

2073 

24 

2074 

25 

2075 

26 

2076 

27 

2077 

28 

2078 

29 

2079 

C 

2080 

2081 

2082 

2083 

2084 

C 

2085 

2086 
2087 

C 

2088 

2089 

2090 

2091 

8 

2092 

7 

2093 

6 

2094 

2095 

2096 

C 

2097 

c - 

2098 

C 

2099 

2  100 

c 

COPY  APPROPRIATE  ARRAY  TO  TEM 1 

I F ( I CODE . EQ . 3 )  CALL  COPYR ( TEMP , 1 , TEM 1 , 1 , NDX*NDY ) 

I F ( ICODE . EQ . 4  )  CALL  COPYR ( POTENT ,  1 , TEM 1 ,  1 , NDX+NDY ) 

IF( ICODE . E0.5)  CALL  COPYR ( QTHERM , 1 , TEM 1 , 1 , NDX*NDY ) 

I F ( ICODE . EO . 6 )  CALL  COPYR ( ELECON , 1 , TEM 1 , 1 , NDX*NDY ) 

DO  INTERPOLATION 

RESULTS  STORED  IN  TEM3 ,  X  COORDINATES  IN  T1,  Y  COORDINATES  IN  T2 
DO  1  J  = 1 , NPY 

Y=YMIN+ (J-1 )*SY 
T2 ( J  )  = Y 
DO  2  1=1 , NPX 

X=XMIN+( 1-1 ) *SX 
T  1  (  I  )  =  X 

TEM3( I , J)=ALINTP(X, Y.TEM1 , XCOORD , YCOORD , NX , NY ) 

CONTINUE 

CONTINUE 

IF  NO  PRINTING  IS  DESIRED,  RETURN 
I F ( IP .NE . 1 )  RETURN 

PRINT  OUT  ARRAY 
PRINT  OUT  HEADING 

I F ( ICODE . EO . 3 )  WRITE( NPRINT , 22 ) 

IF( ICODE .E0.4)  WRITE ( NPRINT , 23 ) 

IF( ICODE . EO. 5)  WR I TE ( NPR I  NT , 24 ) 

I F ( I  CODE . EO . 6 )  WR I TE ( NPR I NT , 25 ) 

FORMAT  STATEMENTS 

FORMAT( 'O' , 10X, 'TEMPERATURE  IN  DEGREES  CELSIUS') 

FORMAT( 'O' , 10X, 'ELECTRICAL  POTENTIAL  IN  VOLTS') 

FORMAT( 'O' , 10X, 'HEATING  RATE  IN  J0ULES/SEC0ND~METER**3 ' ) 
FORMAT( 'O' , 10X, 'ELECTRICAL  CONDUCTIVITY  IN  MHOS/METER') 

FORMAT ( ' 0  COLUMN  ',10(8X,I2)) 

FORMAT ( ' 0  X  (M)' , 10( 3X , F7 . 3 ) ) 

FORMAT ( '  Y  (M)') 

FORMAT ( '  ' .F7.3.5X, 10(2X,F8. 1)) 

DO  6  11=1 .NPX, 10 

I T  =MINO( I  1+9, NPX) 

DO  7  J2= 1 , NPY, 50 

JI =MINO( U2+49 , NPY ) 

PRINT  PAGE  HEADING 

WRITE( NPRINT, 27)  (T  1  ( I  )  , I  =  1 1 , IT ) 

WRITE ( NPRINT , 28 ) 

PRINT  EACH  ROW 

DO  8  JM= J2 , JI 
J  =  NP Y+ 1  - JM 

WRITE (NPRINT, 29)  T2( J) . (TEM3( I, J). 1=11, IT) 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 


SUBROUTINE  OUTR(A.B) 


2101 

C 

3.5 

OUTPUT  A  REAL  VARIABLE  WITH  A  LABEL 

2102 

C 

2103 

c/ 

INSERT  COMBAS 

2104 

REAL  A ( 8 ) , B 

2105 

WRITE ( NPR INT . 20 )  ( A ( I  )  , I  =  1 , 8 ) . B 

2106 

20 

FORMAT ( 4X , 8A4 , '  =',1PE12.4) 

2107 

RETURN 

2108 

END 

2109 

2110 

c 

c _ 

V 

C 

SUBROUTINE  OUTI(A.B) 

2111 

2112 

2113 

C 

2  1  14 

c 

3.6 

OUTPUT  A  LABEL  AND  AN  INTEGER  VARIABLE 

2115 

c 

21  16 

c/ 

INSERT  COMBAS 

2117 

REAL  A ( 8 ) 

21  18 

INTEGER  B 

21  19 

WRITE ( NPR INT , 20 )  ( A ( I  )  , I  =  1 , 8 ) , B 

2120 

20 

FORMAT ( 4X , 8A4 . '  =' . I 12) 

2121 

RETURN 

2122 

2123 

C 

END 

2124 

c _ 

c 

2125 

2126 

2127 

c 

SUBROUTINE  OUTH(A,B) 

2128 

c 

3 . 7 

OUTPUT  A  LABEL  AND  A  HOLLERITH  VARABLE . 

2129 

c 

2130 

c/ 

INSERT  COMBAS 

2131 

REAL  A ( 8 ) 

2132 

REAL  B 

2133 

WRITE (NPR I NT ,20)  ( A (  I  )  ,  I  =  1 . 8  )  . B 

2134 

20 

FORMAT ( 4X , 8A4 , '  = ' , 8X , A4 ) 

2135 

RETURN 

2136 

END 

2137 

2138 

C 

c  — 

2139 

2140 

c 

SUBROUTINE  OUTTAP ( ICODE ) 

2141 

c 

m 

2142 

c 

3 . 8 

OUTPUT  RUN  INFORMATION  TO  MAGNETIC  TAPE 

2143 

c 

2144 

c/ 

INSERT  COMBAS 

2145 

c/ 

INSERT  COMGLO 

2146 

c/ 

INSERT  COMELE 

2147 

c/ 

INSERT  COMTHM 

2148 

c/ 

INSERT  COMCON 

2149 

c/ 

INSERT  COMDIM 

2150 

2151 

c 

GOTO  (1,2)  ,  ICODE 

2152 

c 

2153 

c 

INITIAL  STORAGE 

2154 

1 

WRITE (NONLIN)  ELEGEO , THMGEO , DELTAX . DELTA Y , THICK , 

2  155 

M 

XCOORD , YCOORD , XMIN , XMAX , YMIN , YMAX , 

ELEVOL, 

2156 

H 

ELEALP.ELEBET.EPSELE, THMCON, THMCAP 

, THMTEM , 

2157 

H 

T INI T . NX . NY , NALPHA , MMAX , NREG , NGEO , 

NSPO , 

2  158 

H 

NSTO.NPX.NPY.NIJ.NCI 

2159 

NREC  =NREC+ 1 

2160 

RETURN 

r 


\ 


2161 

2162 

2163 

2164 

2  165 

2166 

2167 

2168 

2169 

2170 

2171 

2172 

2173 

2174 

2175 

2176 

2177 

2178 

2179 

2180 

2181 

2182 

2183 

2  184 

2185 

2186 

2187 

2188 

2189 

2190 

2191 

2192 

2193 

2194 

2195 

2196 

2197 

2198 

2199 

2200 

2201 

2202 

2203 

2204 

2205 

2206 

2207 

2208 

2209 

2210 

221  1 

2212 

2213 

2214 

2215 

2216 

22  17 

2218 

2219 

2220 
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c 

c 

2 


c 

c 

c 


PERIODIC  STORAGE 

WRITE (NONLIN)  NSTEP , TIME , POTENT , QTHERM , TEMP , TEENG , TENG . 
*  TQENG , OENG , NREC , DTI ME 

NREC=NREC+ 1 
RETURN 

END 


SUBROUTINE  TESEND 
C 

C  4.1  TEST  FOR  COMPLETION  OF  RUN 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMGLO 
C/  INSERT  COMCON 
C/  INSERT  COMDIV 
C 

I F ( NFRROR . NE . 0 )  CALL  ENDRUN 
I F ( NSTEP . GE . NRUN )  NLEND  =  .TRUE. 
IF(TIME.GE. (HTIME+CTIME ) )  NLEND  =  .TRUE. 
RETURN 
END 
C 

c - 

c 

SUBROUTINE  ENDRUN 
C 

C  4.2  TERMINATE  THE  RUN 
C 

C/  INSERT  COMBAS 
C/  INSERT  COMDDP 
C/  INSERT  COMDIM 
C 

C  CHECK  IF  A  NORMAL  TERMINATION 

IF ( .NOT .NLEND)  GO  TO  100 
CALL  BLINES(2) 

CALLMESAGE ( 48H  4.2  TERMINATE  THE  RUN 

WRITE ( 18,20) 

20  FORMAT ( ' + ' , 50X , ' ***  NORMAL  EXIT') 

CALL  BLINES ( 2 ) 

CALL  DAYTIM 
CALL  RUNTIM 
CALL  PAGE 
STOP 
C 

c  ERROR  WAS  FOUND  SOMEWHERE,  ABNORMAL  TERMINATION 

100  CONTINUE 

CALL  BLINES ( 3 ) 


CALLMESAGE ( 48H  4.2  ABNORMAL  EXIT  ) 

I F ( NERROR . EO. 1 ) 

*  CALLMESAGE ( 48H  ***  SOLUTION  FAILED  TO  CONVERGE  IN  ELEPOT  ) 

I F ( NERROR . EO . 2 ) 

H  CALLMESAGE ( 48H  ***  GEOMETRIC  ERROR  DETECTED  IN  ROUTINE  OCALC) 

I F ( NERROR . EO . 3 ) 

M  CALLMESAGE ( 48H  ***  ERROR  RETURN  FROM  ROUTINE  TCALC  ) 

I F ( NERROR . EO . 4 ) 

H  CALLMESAGE ( 4 8H  ***  ERROR  RETURN  FROM  SUBROUTINE  ELECOF  ) 


CALL  OUT GRD ( 4 ) 


) 


2221 

2222 

2223 

2224 

2225 

2226 

2227 

2228 

2229 

2230 

2231 

2232 

2233 

2234 

2235 

2236 

2237 

2238 

2239 

2240 

2241 

2242 

2243 

2244 

2245 

2246 

2247 

2248 

2249 

2250 

2251 

2252 

2253 

2254 

2255 

2256 

2257 

2258 

2259 

2260 

2261 

2262 

2263 

2264 

2265 

2266 

2267 

2268 

2269 

2270 

2271 

2272 

2273 

2274 

2275 

2276 

2277 

2278 

2279 

2280 
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CALL  0UTGRD( 5 ) 
CALL  OUTGRD ( 3 ) 

C 

C  PRINT  THE  FINAL  TIMES 

CALL  BLINES( 2 ) 

CALL  DAYTIM 
CALL  RUNTIM 
CALL  PAGE 
C 


C 

c- 

c 


STOP  2230 
END 


SUBROUTINE  MESAGE ( KMESS ) 


C  U.1  PRINT  48-CHARACTER  MESSAGE  ON  OUTPUT  CHANNEL 
C 

C  VERSION  IB  17/12/73  KVR/MHH 

C 

C/  INSERT  COMBAS 

DIMENSION  KMESS ( 12) 

C 

WRITE(NOUT ,9900)  ( KMESS ( J ) . U= 1 ,12) 

C 

RETURN 

9900  F0RMAT(4X, 12A4) 

END 


C 

C 

C 


SUBROUTINE  PAGE 
C 

C  U.2  FETCH  NEW  PAGE  ON  OUTPUT  CHANNEL 
C 

C  VERSION  IB  17/12/73 

C 

C/  INSERT  COMBAS 

WRITE ( NOUT , 9900 ) 

C 


RETURN 

9900  FORMAT ( 1H 1 ) 

END 


KVR/MHH 


C 

c - 

c 

SUBROUTINE  BLINES(K) 


C 

C  U.3  INSERT  BLANK  LINES  ON  OUTPUT 
C 

C  VERSION  IB  17/12/73 

C 

C/  INSERT  COMBAS 
C 


100 

DO  100  J=  1  , K 

WR I T  E ( NOUT , 9900) 

9900 

RETURN 

FORMAT ( 1H  ) 

END 

CHANNEL 

KVR/MHH 


CULHAM 


CULHAM 


CULHAM 
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SUBROUTINE  RUNTIM 

U . 1 2  UPDATE  CPU  TIME  (MINS)  AND  PRINT  IT 

VERSION  2B  17/12/73  KVR/MHH  '  CULHAM 

/  INSERT  COMBAS 

TIME ( 1 , 0 , ICPUT )  IS  A  MTS  LIBRARY  FUNCTION  GIVING  THE  CPU-TIME 
USED  SO  FAR  (ICPUT  IN  MSECS) 

CALL  TIME ( 1 ,0, ICPUT) 

I  SEC  =  ICPUT/ 1000 
MIN  =  I  SEC/60 
I  SEC  =  I  SEC  -  60*MI N 
SEC  =  I  SEC 

CPTIME  =  MIN  +  SEC/60. 

WRITE ( NOUT , 9900 )  MIN.ISEC 

RETURN 

9900  FORMAT ( 5X , 22HCPU  TIME  USED  SO  FAR  =,I4,5H  MINS, 14, 5H  SECS) 
END 


SUBROUTINE  DAYTIM 


U . 1 3  PRINT  DATE  AND  TIME 


VERSION  2B 

17/12/73 

KVR/MHH 

CULHAM 

INSERT  COMBAS 

TIME  IS  AN  MTS 
DIMENSION 

LIBRARY  ROUTINE 
DATE (4),TIIME(2) 

WHICH  RETURNS  THE 

DATE  AND  TIME 

CALL  TIME ( 2 1 , 0 , DATE ) 

CALL  TIME(4,0,TIIME) 

WRITE(NOUT , 9900)  DATE , TI IME 
RETURN 

9900  F0RMAT(5X,4A4,5X,2A4) 

END 


C 

C 

C 

C 

C 

C 


C 


SUBROUTINE  RESETR( PA , KDIM , PVALUE ) 

U.14  RESET  REAL  ARRAY  TO  SPECIFIED  VALUE 

VERSION  IB  17/12/73  KVR/MHH 

DIMENSION  PA(KDIM) 

DO  100  U= 1 , KDIM 
PA ( J )  =  PVALUE 
100  CONTINUE 

RETURN 

END 


CULHAM 


X 


sJ 
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SUBROUTINE  RESETI ( KA , KDIM , K VALUE ) 

U . 1 5  RESET  INTEGER  ARRAY  TO  SPECIFIED  VALUE 

VERSION  IB  17/12/73  KVR/MHH  CULHAM 

DIMENSION  KA ( KDIM) 

DO  100  J=1 , KDIM 
KA ( J )  =  KVALUE 
100  CONTINUE 

RETURN 

END 


SUBROUTINE  RESETH( KA , KDIM , KVALUE ) 

U . 1 6  RESET  HOLLERITH  ARRAY  TO  SPECIFIED  VALUE 

VERSION  IB  17/12/73  KVR/MHH  CULHAM 

DIMENSION  KA(KDIM),  KVALUE(I) 

DO  100  J  = 1 .KDIM 
KA ( J )  =  KVALUE ( 1 ) 

100  CONTINUE 

RETURN 

END 


SUBROUTINE  JOBTIM( PTIME ) 


U . 1 7  FETCH  ALLOCATED  JOBTIME  (MINS) 

VERSION  2B  17/12/73  KVR/MHH  CULHAM 


SET  CPU  TIME  COUNTER  TO  0. 

CALL  TIME(0.0, JUNK) 

GUINFO( 78 , INF0 1 )  IS  AN  MTS  LIBRARY  FUNCTION  WHICH  GIVES  THE  GLOBAL 
TIME  REMAINING 

GUINFO( 86.INF02)  IS  AN  MTS  LIBRARY  FUNCTION  WHICH  GIVES  THE  LOCAL 
TIME  REMAINING 

INF01  AND  INF02  ARE  IN  CPU  TIME  UNITS.  ONE  TIME  UNIT 
EQUALS  13.0208333333  MICROSECS.  CPTIME  IS  IN  MINS. 

CALL  GUINFO( 78.INF01) 

CALL  GUINFO( 86, INF02) 

INFO  =  MINO( INF01 , INF02) 

PTIME  =  2 . 1701 389* INFO 


RETURN 

END 


'■fi  wiTt m 
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SUBROUTINE  RESETL ( KLA , KDIM , KLVAL ) 

U . 19  RESET  LOGICAL  ARRAY  TO  SPECIFIED  VALUE 

VERSION  IB  17/12/73  KVR/MHH 

LOGICAL  KLA,  KLVAL 

DIMENSION  KLA ( KDIM ) 

DO  100  0=1 , KDIM 
KLA ( 0 )  =  KLVAL 
100  CONTINUE 

RETURN 

END 


SUBROUTINE  SCALER ( PA , KDIM , PC ) 

U . 2 1  SCALE  A  REAL  ARRAY  BY  A  REAL  VALUE 

VERSION  IB  17/12/73  KVR/MHH 

DIMENSION  PA ( KDIM ) 

DO  100  0=1 .KDIM 
PA ( 0 )  =  PA ( 0 )  *  PC 
100  CONTINUE 

RETURN 

END 


SUBROUTINE  C0PYR(PA1 , K 1 , PA2 , K2 , KDIM ) 

U . 23  COPY  ONE  REAL  ARRAY  INTO  ANOTHER 

VERSION  IB  17/12/73  KVR/MHH 

DIMENSION  PA  1 ( KDIM ) , PA2 ( KDIM ) 


DO  100  0= 1 , KDIM 
I  1  =  K1  +  0  -  1 
12  =  K2  +  0  -  1 
PA2 (12)  =  PA  1(11) 
100  CONTINUE 

RETURN 

END 


MODULE  COMBAS 


1 . 1  BASIC 

COMMON  /COMBAS/ 

R  ALTIME , 

I  LABEL  1  , 

I  LABEL5 , 

I  NDIARY. 


SYSTEM  PARAMETERS 


CPTIME 

LABEL2 

LABEL6 

NIN 


STIME , 
LABEL3 , 
LABEL7 , 
NLEDGE . 


LABEL4  . 
LABEL8  , 
NONLIN, 


CULHAM 


CULHAM 


CULHAM 


9 

10 

1  1 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

4  1 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 


1  1  1 


I 

I 

L 


NOUT ,  NPRINT,  NPUNCH,  NRE AD , 

NREC,  NRESUM,  NRUN,  NSTEP, 

NLEND,  NLRES 


LOGICAL  NLEND, NLRES 
C 


c/ 

c 

c 


c 


c 


INTEGERS  LABEL  1(  12)  ,LABEL2(  12)  ,LABEL3(  12)  ,LABEL4(  12)  , 
I  LABEL5( 12) . LABEL6 ( 12) , LABEL7 ( 12) ,LABEL8( 12) 

MODULE  COMDDP 


1.9  DEVELOPMENT  AND 

COMMON  /COMDDP/ 

I  MAXDUM ,  MXDUMP , 

I  NPDUMP ,  NPOINT , 

L  NLCHED ,  NLHEAD , 

L  NL0MT3 ,  NLREPT 


DIAGNOSTIC  PARAMETERS 

NADUMP,  NCLASS , 

NSUB ,  NVDUMP , 

NL0MT1,  NL0MT2 , 


LOGICAL  NLCHED, NLHEAD(9) , NLOMT 1 (50) ,NL0MT2(50) ,NL0MT3(50) , 
L  NLREPT 


C/ 

C 

C 


C 


C/ 

C 

C 


C 


C/ 

C 

C 


C 


c/ 

c 

c 


INTEGERS  NADUMP  (  20 )  ,NPDUMP(20)  ,NVDUMP(20) 

MODULE  COMGLO 

1 . 2  GLOBAL  VARIABLES 

COMMON  /COMGLO/ 

P  ELECON, OTHERM, POTENT , TEMP , 

P  DELTAX.DELTAY .THICK, DTI ME .TIME , 

P  XCOORD . YCOORD . XMIN , XMAX , YMIN , YMAX , 

P  ELEGEO , THMGEO 

REAL  ELECON( 50 , 50) , OTHERM ( 50 . 50 ) ,DELTAX( 50 ) . 

P  POTENT ( 50 , 50) , TEMP ( 50 , 50 ) ,DELTAY( 50) , 

P  XCOORD ( 50) , YCOORD ( 50 ) 

INTEGERM  ELEGE0(50, 50)  ,  THMGEO(  50, 50) 

MODULE  COMELE 

1.3  ELECTRICAL  CONDUCTIVITY  AND  DIFFERENCE  EOU.  COEFFICENTS 

COMMON  /COMELE/ 

P  ECXP, ECXM, ECYP, ECYM, ERHS. EXMXP, EYMYP, 

P  ELEVOL, ELEALP, ELEBET, ELEPAR, EPSELE 

REAL  ECXP ( 50 , 50 ) , ECXM( 50 , 50 ) , ECYP ( 50 , 50 ) , ECYM( 50 , 50 ) , 
P  ELEVOL ( 10) , ELEALP( 10) . ELEBET ( 10) , ELEPAR ( 50) , 

P  ERHS ( 50 , 50 ) , EXMXP ( 50 ,50),EYMYP( 50 , 50 ) 

MODULE  COMTHM 

1.4  THERMAL  CONDUCTIVITY  AND  HEAT  CAPACITY 

COMMON  /COMTHM/ 

p  THMCON , THMCAP , THMTEM , TCXP . TCXM , TCYP , TCYM . 

P  TCX , TCY , TCRHS 

REAL  THMCON ( 10) ,THMCAP( 10) ,THMTEM( 10) 

REAL  TCXP ( 50 , 50 ),TCXM(50,50) , TCYP( 50. 50) . TCYM( 50, 50) 
REAL  TCX ( 50 , 50 ) , TCY ( 50 . 50 ) , TCRHS ( 50 , 50 ) 

MODULE  COMTEM 

1 . 5  TEMPERARY  ARRAYS 

COMMON  /COMTEM/ 

M  TEM1 .TEM2.TEM3.T1 .T2.T3.T4.T5 


C 


r  ^ 


€9 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

OF 
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REAL  TEM1(50,50) ,TEM2( 50,50 ) ,TEM3( 50,50) ,T1( 50) , 

#  T2 ( 50 ) , T3 ( 50 ) , T4 ( 50 ) , T5 ( 50 ) 

C/  MODULE  COMDIM 

C 

C  1.6  DIMENSION  AND  COUNTER  VARIABLES 
COMMON  /COMDIM/ 

#  NX , NY , NEOU , NALPHA .MITER . MMAX . NREG . NGEO . NERROR , NTYPE . 

#  NSPO , NPX , NPY , NI J , NCI , NDX , NDY , NSTO 
C 

INTEGERS  NX  .  NY  ,  NDX  .  NDY  .  NEOU  ,  NALPHA  ,  MITER  .  MMAX 
INTEGERM  NREG ,  NGEO  .  NERROR  .  NTYPE  ,  NSPO  .  NPX  ,  NPY 
INTEGER*4  NIJ, NCI. NSTO 
C/  MODULE  COMCON 
C 

C  1.7  CONTROL  THE  HEATING 
COMMON  /COMCON/ 

H  HTIME , CTIME . CCUR . CPOW , DTEMP . DELT . CUR . 

#  VOLTS. RESI ST. POWER. CVOL , TINIT , DEENG , TEENG , 

M  DOENG , OENG , TOENG , TTENG , TENG , SCALE 

C 

REAL  OENG( 10) ,TENG( 10) 

FILE 


t 


X 


1 
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10 

1  1 
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21 
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23 
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26 

27 

28 

29 
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31 

32 

33 
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Appendix  II  -  Sample  Output  from  MEGAERA 


&INMOD  NRUN  = 100  SEND 
&LABELS 

LABEL1='MEGAERA  RUN  NO.  35  PARALLEL  ELECTRODES', 

LABEL2= 'OVERBURDEN  CONDUCTIVITY  RATIO  OS/OB  =  1.0', 
LABEL4= ' GEOMETRY  SET  UP  AND  ENTERED  BY  ALLAN  HI  EBERT ' , 

&END 

SINPUT1  NX=50 , NY=50 , 

DELTAX=50* 1 .5,DELTAY=8*5.0,2.5,3*2.0, 1.5,25*1.2, 

1.5, 3*2. 0,2. 5, 6*5.0, 10.0,NREG=5, 

SEND 

SREGION  MINI  =  1 , MAX  I =50 , MINJ= 1 4 , MAXd  =  38 , 

ETYPE= 'DOMA ' , E ALPHA =2 . 29E-2 , EBETA= 1 .E-3, 

TTYPE= ' COND ' ,THMRC=1 .6E6,THMK=1 .8, 

SEND 

&REGION  MINI = 1 , MAXI =50, MINJ= 1 ,MAXd=13, 

ETYPE  = ' DOMA ' , E ALPHA  =  2 . 29E-2 , EBETA= 1 .OE-3, 

TTYPE= ' COND ' , THMRC= 1 . 6E6 , THMK= 1.8, 

&END 

SREGION  MINI = 1 , MAXI =50 , MINJ=39 , MAXJ=50 , 

ETYPE= ' DOMA ' , E ALPHA =2 .  29E~2,EBETA=1 .OE-3, 

TTYPE= ' COND ' ,THMRC=1 ,6E6,THMK=1 .8, 

&END 

SREGION  MINI = 1 ,MAXI=1 , MINd=24 , MAXU=29 , 

ETYPE  = ' ELEC ' ,V0LTS=1OOO.  , 

TTYPE  = ' COND ' ,THMRC=1 . E6 , THMK  =  3 . O, 

SEND 

SREGION  MINI =50 , MAXI =50 , MINJ=24 , MAXU=29 , 

ETYPE= ' ELEC ' ,V0LTS=0. , 

TTYPE =' COND' ,THMRC=1 . E6 , THMK=3 . O , 

SEND 

SINPUT2 

NSPO= 10,XMIN=1 . 0 , XMAX=49 . 0 , NPX= 1 7  , 

YMIN=40 , YMAX=85 , NPY= 10 , 

SEND 

SINPUT3 

TINIT=15.  ,DTIME=1 .E5,HTIME  =  3.  1 54E7 , NTYPE  =  3 , 

CPOW= 12 . E3 , CVOL= 1000 . , DTEMP= . 05 , NCI =20 , 

MMAX=350, EPSELE= .02 . 

SEND 

SINPLOT 

LABEL1  =  'PARALLEL  PLATE  RUN  NO.  35  80/ 1 1 /25 ' , NLAB 1 =38  , 

LABEL2= 'CONDUCTIVITY  RATIO  OB/OS  =  1.0/  1 . O ' . NLAB2=38 , 
NCPLTS= 1 ,CPTIME  =  3 . 15E7 , NCC0DE  =  3  , 

CXMI N= 1 . 0, CXMAX=74 . O, CYMIN=30 . ,CYMAX=100. , 

NPPLTS= 1 ,PPTIME=3. 1 5E7 , NPCODE =7 , 

PXMIN= 1 . 0 , PXMAX=99 .0,PYMIN=15. ,PYMAX=1 15. , 

SEND 


MEGAERA  RUN  NO.  35  PARALLEL  ELECTRODES 
OVERBURDEN  CONDUCTIVITY  RATIO  OS/OB  »  1 
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